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ABSTRACT 


This synopsis of compressible fluid dynamics was read in a 
shortened form before the Seventh International Congress of 
Applied Mechanics. It treats successively, with full references, 
the foundations of the subject, the boundary layer, isentropic 
flow, plane waves, and the differences between steady sub- and 
supersonic flow; the perturbation methods of solution, including 
the linearized theory of supersonic flow in its many ramifications, 
are described and compared with more exact work in certain 
cases, and the theory of characteristics is outlined. The memo- 
randum ends with a full discussion of the uses of the hodograph 


transformation. 


INTRODUCTION 


BEHAVIOR OF GASES in high-speed flow is of in- 
creasing engineering importance. Experiments on 
it are too often costly failures and the cost may include 
human lives; hence, methods for predicting as much as 
possible in advance are extremely desirable, and it is 
these I shall attempt to describe in this lecture. 

The foundation of our approach is the Kinetic Theory 
of Gases, which explains their behavior in terms of the 
“tyranny of large numbers.’’ So well are its more care- 
ful predictions in agreement with experiments, where 
any have been made, that we are content to accept 
them also when no such checks are available. By its 
means a simplified treatment of gaseous flow is adopted, 
which assigns to each point of the gas at each time a 
density, temperature, velocity, and stress tensor, with a 
static pressure p which is the mean of the normal 
stresses in three mutually perpendicular directions. 
Each of these quantities is a continuous function in 
space and time, being defined as a local average of an 
appropriate physical property of the molecules. By an 
averaging of the laws of motion of the molecules, basic 
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postulates for gas dynamics are obtained. There is a 
relation®* between the rate of strain and the stress tensor 
involving a single parameter yu, the viscosity, varying 
from point to point but depending for a given gas on the 
temperature in a manner usually approximated ade- 
quately by Sutherland’s formula? p«7°*/(T + B), 
where T is the absolute temperature and B = 114° for 
air. There is an equation of state connecting density, 
temperature, and pressure; an equation of conservation 
of mass; and a law of transport of temperature*® con- 
taining a parameter k, the conductivity, also dependent 
on temperature being usually closely proportional to u. 
For most gases w and & are small enough to render 
negligible the tangential stresses, the differences be- 
tween the normal stresses and the heat transfer, unless 
either the velocity gradient or the temperature gradient 
is large; but these two conditions are practically 
equivalent, and the effects of viscosity and conductivity 
are felt together and are usually of comparable magni- 
tude. 

The conditions in question will arise in the first in- 
stance near any solid boundary, where the postulate of 
continuous velocity implies a rapid drop to zero 
through some ‘‘boundary layer.’’ They will also arise 
in a region where compression waves are piling up on 
one another either in space or time when a régime is set 
up called a shock wave with a large change in velocity, 
temperature, density, and pressure across a surface of 
microscopic thickness, with considerable viscous dissi- 
pation of energy and conduction of heat within it. The 
gas gains entropy in passing through the shock wave so 
that less energy is available behind it for appearance in a 
kinetic form. 

The isentropic flow when both yu and & are identically 
zero is a reasonable approximation to the truth in the 
initial stages of any disturbance from a uniform state 
and also later, provided that the parts of the gas whose 
energy is dissipated have a fairly narrow outlet to in- 
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finity. Thus boundary-layer separation must be 
avoided; we know that at low speeds this is so if any 
retardation of the main stream at the surface is not too 
great. Investigations extending this conclusion to 
higher speeds are in their infancy, but at any rate ex- 
periments indicate that dissipation phenomena are 
negligible in sufficiently many problems to make theo- 
retical work on this hypothesis useful. Even more 
problems can be studied if all shock waves be replaced 
by surfaces of zero thickness across which absolute dis- 
continuities in physical quantities occur governed by 
relations due to Rankine* ?® and rediscovered by 
Hugoniot. When this is done, the energy, though con- 
stant behind the shock wave for each fluid element, may 
differ for different elements. 

Among investigations not making these simplifica- 
tions is the important work on the shock wave of finite 
thickness by Rayleigh’! and by Taylor’? using the 
equations of viscous conducting flow. This was carried 
on by Becker'® and Thomas,*” but the thickness de- 
duced was hardly larger than the mean free path of 
molecules in the gas, so a more fundamental approach 
was adopted by Bethe and Teller;* all this work puts 
Rankine’s equations on a firm basis. Again, the lam- 
inar boundary layer in a gas has been investigated by 
many authors. Most have studied it in the case of a 
uniform main stream; Emmons and Brainerd® used 
the differential analyzer on this problem, and L. 
Crocco* gave fairly simple solutions either when the 
Prandtl Number is unity (following Busemann,® von 
K4rman and Tsien*) or when the viscosity varies di- 
rectly as the absolute temperature, and he gave a more 
complicated procedure for the general case, concluding 
that the second assumption gives a better approxima- 
tion for air than the first. Recently, Hartree and 
Cope” have applied electronic computing methods to 
the boundary layer in a retarded flow, and Howarth” 
has shown how, by a transformation, the equations 
with any main stream pressure distributions when both 
of Crocco’s simplifying assumptions are made do not 
become much more complicated than in the low-speed 
case. He solves various problems and draws general 
conclusions that, if a laminar boundary layer of com- 
pressible fluid be compared with an incompressible one 
with the same main stream velocity distribution and 
kinematic viscosity, the compressible one has a rather 
larger thickness, an earlier separation, and a reduced 
skin friction coefficient. At Mach Number 1 the 
thickening is about 10 per cent; at Mach Number 2, 
about 50 per cent; separation occurs earlier by 30 per 
cent at M = 3.16, occurring, in fact, for a linear adverse 
gradient when has dropped to 2.4. Other investiga- 
tions admitting dissipative phenomena have been 
mostly semiempirical and ad hoc, such as the postulates 
on shock wave-boundary layer interaction, the ingen- 
ious method of Cope’ for estimating the base drag of a 
shell, and the theory of condensation shocks, which 
shows, as is observed, that condensation of water vapor 


with its resulting gain of heat may happen in a thin 
region analogous to a shock wave but with velocity 
changes opposite in sign. 


ISENTROPIC FLOW 


The rest of what I shall say concerns adiabatic flow— 
that is, flow without heat transfer between fluid ele- 
ments. It may be divided into the isentropic flow, 
where shock waves are absent or uniform in character, 
and the anisentropic flow behind a nonuniform shock 
wave. Isentropic flow implies a functional relationship 
between pressure and density; for a perfect gas with 
constant specific heats, this is p/p’ = constant, where 4 
is the ratio of the specific heats. For other gases it is 
obtained from the differential equation 


dp/dv = —y([(OT Ov) »/ (07 /dp),] 


where v = 1/p. Though much I shall say is true for any 
isentropic motion, few calculations have been attempted 
except on the assumption p/p” = constant. The 
assumption is reasonable for air, but the value of y to be 
adopted is contested, for it is difficult to measure. A 
committee of the British A.R.C. has adopted a resolu- 
tion that the value y = 1.400 should be used as much as 
possible since experiments indicate that y is nearer to 
this value, especially in the conditions of high-speed 
flow applicable in aeronautics, than to the other favored 
by investigators, which is y = 1.405. 

For isentropic gas flow Kelvin’s result’ *? holds that 
the circulation in a circuit moving with the fluid re- 
mains the same as long as no viscous action takes place 
at points of the circuit. There is no restriction on what 
may happen inside the circuit, but it must cross no dis- 
continuities of velocity such as vortex sheets or shock 
waves if the theorem is to be true; when it does cross a 
shock wave, the rate of change of circulation is not 
necessarily only of the order of the variations of entropy 
in the field (which may be extremely small) but can be 
of the order of the velocity change across the shock 
wave. From Kelvin’s theorem are deduced Helm- 
holtz’s result,®: 7? that vortex lines move with the fluid, 
and Prandtl’s arguments’ establishing the general form 
of the trailing vortex theory of flow past a wing, which 
therefore hold rigorously below the speed at which a 
shock wave first appears. Another consequence (stated 
first by Lagrange’) is that portions of the fluid moving 
irrotationally continue to do so unless they are sub- 
jected to the action of viscosity or conductivity. 
Again, Joukowski’s empirical law for airfoils,* that the 
circulation is such as will make the flow come off 
smoothly at the trailing edge, appears to be true with 
much the same qualifications as for low-speed flow; 
while his theorem, that in two-dimensional flow the lift 
is pU times the circulation, is true in isentropic irrota- 
tional flow if p the density as well as U the velocity is 
given its main stream value. But again, when shock 
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waves are present, the error may be of greater order 
than the variations of entropy in the field. 

In the irrotational motion that Lagrange’s result 
assures for us under certain conditions, a velocity 
potential ¢ exists, and the momentum equation can be 
integrated to yield Bernoulli’s equation” 


Pay 2d 1 
~ re [% += 9° = constant (1) 


< 


where a is the local velocity of sound and a? = dp/dp. 
[When p <p”, the integral can be written a*/(y — 1).] 
can be chosen so that the constant on the right does 
not vary with time. The equation of conservation of 


mass?’ is 


Op/Ot + div(p grad ¢) = 0 


or 
LO l 
V%>=- ( ~4 + - grad p grad ‘) 
pot p 
But 
dp/p = —a~*d[(O0o/Ot) + (1/2)q?] 
so that 


LV °h = 0°p/Ol? + 0(1/2q") /Ot + 
grad @ grad [(0¢/Odt) + (1/2)q?] 


ome. oa « ad - ca. - 
av -¢ = wv + OE + grad ¢ grac (50) (2) 
Since qg? is (grad ¢)? and a? is expressible in terms of 
0¢/Of and q?, this is a single partial differential equation 
for @ in four variables. When a is large compared both 
with the fluid velocities and with the speed at which the 
field is changing, the equation is approximable as V *¢ = 
0. I need say little about this case except to remark 
that the motion then depends only on the instantaneous 
motion of the boundary, not on its previous movements, 
When rapid changes in the field are occurring, the first 
term 0°¢/0/? on the right must be taken into account, 
and superimposed on the previous flow we may have 
sound waves moving about. Finally, when large 
Mach Numbers g/a appear, the whole motion of the gas 
becomes dependent on the life history of the surfaces 
producing it, and ¢ satisfies a nonlinear equation. 
However, when uniform boundary conditions have 
obtained for some time, a steady flow with @ inde- 
pendent of ¢ may result. Other flows including those 
with sound waves and turbulence are called unsteady, 
but I am relieved from discussing turbulence here and 
the theory of sound is not high-speed gas flow. Work- 
ers on laminar unsteady high-speed flow have studied 
problems where there is one independent variable in- 
stead of four, as in Taylor’s study® of the spherical 
shock wave, or two variables as in the problem of rec- 
tiiinear motion I shall now discuss, or they have 


adopted perturbation methods which will be mentioned 
later. 


ONE-DIMENSIONAL FLOW 


The equation of one-dimensional isentropic flow is 
homogeneous in the three second derivatives of 9. 
Such an equation becomes linear” if the Legendre 
potential 


w = x0o/dx + 106/dt — 
be used as independent variable and expressed in terms 


of 0¢/Ox and 0¢/dt. It becomes simpler if the variable 
b = 2a/(y — 1) be used instead of 0¢/02, when we have 


db? ' y—1b2b~ dq? 





Ow 3-—yvyldw dw ; 
to. (3) 


with x, ¢ given by 


t = —2(y — 1)—'b-'(0w/0dD) 
x = gt + (Ow/dg) 


Even this equation is complicated in general, and the 
best approach is by the method of characteristics, which 
Riemann invented‘ for this specific problem. This 
approach, described later on, is also applicable to the 
analogous problems with cylindrical or spherical sym- 
metry. The plane waves problem simplifies, however, 
for a progressive wave in one direction only, when the 
general solution is 


a = a) + (1/2)(7 — Deg, q=f(x-(a+ qt] (4 


This is substantially due to Earnshaw;;’ it is easily veri- 
fied to be a solution, and the function f can be found 
from initial conditions. From it is readily deduced"! * 
that in a compression wave the positions where two 
given values of the velocity occur become closer and 
closer, so that the wave becomes more intense, until 
finally the distance between the positions drops through 
zero and becomes negative, an impossible configuration, 
and it must be assumed that a shock wave has super- 
seded the previous isentropic flow. Eq. (3) also simpli- 


fies when 
(3 — y)/(y — 1) = 2m 


an even integer, so that y = 3, 5/3, 7/5, 9/7, etc. The 
first gives the ordinary wave equation for w but corre- 
sponds to no known gas. The second gives the equa- 
tion of spherical wave motion with its simple general 
solution; this value holds for monatomic gases and was 
used by Martin® in his study of the interaction of waves 
of expansion with shock waves. For m 2 1 the general 
solution is 


(02/dg? — d2/db%)"—"{b—[f(q + 6) + (gq — d)]} 


The value, y = 11/9, m = 4, was used by Love and 
Pidduck" in their work on the reflection of waves in gun 


barrels. In practice, the variables, r = (1/2)(q + 4), 
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s = (1/2)(q — b), are usually most convenient; they 
were introduced by Riemann.‘ 

Steady flow may be obtained in a tunnel or in the 
uniform motion of a body in the air if the body is taken 
as the frame of reference. The role played by the 
critical velocity where the Mach Number is 1 in both 
these configurations is now general knowledge. That 
mass flow per unit area is a maximum when M = 1, by 
Bernoulli’s equation, implies that higher speeds can be 
obtained in a pipe only if it converges and then di- 
verges” and that otherwise an excessive pressure differ- 
ence will produce choking. Again, in the subsonic flow 
past a body, all parts of the fluid are affected to some 
extent by the body’s presence—in fact, if the body hasa 
plane of symmetry perpendicular to the stream, any 
isentropic flow pattern will also have this plane of sym- 
metry. But in supersonic flow, only parts of the fluid 
at any instant are “within earshot’ of the body, and 
these lie in a cone that is more or less rounded at its 
vertex. The body pushes a compression wave in front 
of it, which in the steady state has formed into a shock 
wave. 

The transition to the supersonic régime is begun at a 
subsonic speed by the formation of at least one small 
shock wave in the supersonic region near the surface; 
this shock wave grows and moves toward the front as 
the speed is increased, finally coalescing in most cases 
with another shock wave, which appeared ahead of the 
body when ./ reached 1 and has been moving back as A/ 
increased further. When the first shock wave occurs, 
the total drag can no longer be correlated solely witb 
deficiency of momentum in the wake because of viscous 
action in the boundary layer and of trailing vortices in 
the case of a wing, but it must include also that due to 
the gain in entropy at the shock wave. When J ex- 
ceeds 1, a fourth cause of drag appears, the transport of 
deficiency of momentum to infinity not directly but 
diagonally behind the body in the region immediately 
behind the bow shock wave; however, at this speed the 
induced drag is greatly lessened. Thus, in supersonic 
flow a wake traverse, unless it explored the whole 
Mach cone, could not be used to deduce the total drag. 
Similarly, the lift is not associated primarily with down- 

yash in the wake but with downwash at large distances 
in a diagonal direction behind the bow shock wave. 


PERTURBATION METHODS 


Passing to perturbation methods of solving the non- 
linear equation of steady flow, I begin with that due to 
Janzen'’ and Rayleigh,'* which is useful for calculating 
subsonic flows past thick obstacles. It compares the 
gas flow in question with the flow at the same speed past 
the same obstacle of a liquid or of other gases with 
different velocities of sound. For the liquid the Mach 
Number is zero, and for the other gases the free-stream 
Mach Number J/ takes different values. The potential 
¢ exists for each M and so can reasonably be expanded in 


powers of Mas 
& = do + OM? + gM + ... 


Then ¢ is a harmonic function representing the liquid 
flow past the body, and ¢i, ¢2, etc., satisfy equations of 
Poisson’s type growing in complexity and depending on 
the previous terms in the series. Many of them have 
been calculated for different problems, some of which 
might be useful to experimenters, since in a completely 
subsonic flow the use of ¢ and ¢; alone gives a reason- 
able approximation as Eser shows,** who gives a full 
bibliography. General formulas for determining ¢» have 
been given by Kaplan“ and by Goldstein and myself.*4 
Taylor and Sharman?!*? adopted an electrical method 
of solving Poisson’s equation, working what is effec- 
tively the Janzen-Rayleigh method by successive 
approximations instead of by expansions in series. In 
this method any two-dimensional body shape can be 
used with equal ease, and no step is more laborious than 
the one before. When the boundary was a circle, the 
process converged at JJ = 0.4 but diverged slowly at 
M = 0.5. It is not certain whether it begins to diverge 
at the first appearance of a local supersonic speed in the 
flow; or at the point where isentropic flow becomes im- 
possible; or, as recent evidence indicates, somewhere in 


between. 


LINEARIZED THEORY OF TWO0O- DIMENSIONAL SUBSONIC 
FLow 


More usefully, the flow past a surface with given 
velocity may be compared with the flow past another 
surface similar in all respects except that it is reduced in 
scale by a factor 6 in the direction perpendicular to the 
flow. When 6 = O, the flow is undisturbed with 
velocity, say U, parallel to the z-axis; so for other 6 an 
expansion 


@ = Us + 66, + 6%b2 + 


may be expected. The first two terms will give the 
flow past a thin surface. Substituting this potential in 
the equation of steady flow and arranging both sides in 
powers of 6, the first term is 


9 9 + re) _. Ob 
ao"6V "gh, = U o( ve 21) 


or 

Vo, — 117(07¢,/027) = 0 (5) 
which is known as the linearized equation of isentropic 
gas flow. It can be obtained simply from the basic 
equations by neglecting squares of small quantities 
throughout, but some approach like the foregoing is 
necessary if comparison with more exact solutions is to 
be possible. The equations for ¢2, 3, etc., have the 
same left-hand side as Eq. (5), but the right-hand side 
becomes an increasingly more complicated expression in 
the ¢’s already determined. Boundary conditions at 
the surface must also be expanded in powers of 6. 
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Prandtl*? and Glauert” first used the linearized equa- 
tion in subsonic flow; Goldstein and Young* listed the 
more important results, correctly except for some 
points concerning axisymmetrical flow which Sears® 
has cleared up. In the two-dimensional flow past a 
body (with Joukowski’s condition holding at the trailing 
edge if there is one), the excess pressure at any point of 
the surface is greater than its value in the same flow at 
zero Mach Number by a factor 1 V1 = M?; the lift 
and moment are multiplied by the same factor. It 
tends to infinity as 1J —~> 1 but the lift does not really 
become infinite, for at the same time the value of the 
approximation tends to zero. This is proved by con- 
sidering later terms in the power series in 6. Experi- 
ments show that the lift rises at first but falls again 
when ./ approaches 1. Results also exist comparing 
the flow with that for 1 = 0 past a body of arbitrarily 
reduced thickness. But for flow past a thin body of 
revolution, the comparison, to be correct, must be with 
a body of thickness reduced by a specific factor, 
vi - AM*, and the excess pressure is then 1/(1 — M?) 
times as much. This is essentially because the expan- 
sion of the surface pressure in powers of 6 does not con- 
tain terms simply in 6°, 64, etc., as might be expected for 
this problem, but also contains terms in 6° log 6, etc. 
There is no lift on a streamlined body of revolution at 
yaw, but there is a couple tending to increase the yaw, 
equal to 2p? times the volume of the body times the 
angle of yaw and independent of Mach Number. As 
M —> 1, this does not tend to infinity, but it becomes 
valueless exactly as before. In fact, for JJ > 1, if the 
body is pointed at both ends, the corresponding moment 
is exactly half the above result, and clearly the transi- 
tion will not take place discontinuously. But if it is 
rapid, results of this kind must be taken into account in 
flight through the speed of sound. 


LINEARIZED THEORY OF Two-DIMENSIONAL 
SUPERSONIC FLOW 


Remarking only that the technique for finding further 
terms in the power series in 6 has been expounded by 
Kaplan,** I pass on to the application of the method in 
supersonic flow, which is much easier because, where in 
equations like Laplace’s all the boundary conditions 
affect the solution at each point, equations like the 
wave equation have the property that only conditions 
ahead of the point do. We write the linearized equa- 
tion 

0° 0° 0° 


—- t? 
aes PG OOS 








(6) 


where cot’ = J/* — 1, so that yu is the semiangle of the 
Mach cone in which the disturbance is confined and is 
known as the Mach angle. For two-dimensional 
problems the general solution is 


¢ = f(z — ycot wu) + g(z + y cot u) 


The f term represents a disturbance to the stream prop- 
agating itself along lines z — y cot u = constant and the 
g term along lines z + y cot uw = constant. Both are 
easily obtained in different problems. For an airfoil 
with upper surface y = y,(s) and lower surface y = 
y2(z), the upper boundary condition is 


(0¢/Oy),-, = Uyi'(s) 


and disturbances must be propagated in such a way 
that y increases with z. Hence, ¢ = f(z — y cot yw), 
where 

—cot uf’(z) = Uy;'(z) 


The pressure is given by Bernoulli's equation / (dp/p) 
+ (1/2)g? = constant, which can be approximated as 


[(p = Po) ‘po] + (U0 '0z) = 0 


for most purposes. Thus, the excess pressure on the 
upper surface is pU? tan wy;'(z); similarly, on the lower 
surface, it is —pU* tan wye’(z), and these integrate up to 
give a lift of 2oU? (tan u)h, where h/ is the distance from 
leading edge to trailing edge perpendicular to the 
stream. If a@ is the incidence measured from the no- 
lift angle, this means that C, [lift per unit wing area 
divided by (1/2)pU?] is 4a tan yz, or 4a/V M? — bs 
compared with 2ra/V1 — M?* for thin airfoils at sub- 
sonic speeds. (Neither is valid near JJ = 1.) The 
drag is 


l 
pU* tan wf [yi’2(z) + ye’*(z)]d2 
0 


Its coefficient C, differs from its no-lift value Cp, by 
4o?/V M2 — 1, and Cp, is an integral that, if any 
points on the airfoil are taken as fixed, is a minimum 
when these are joined by straight lines. This suggests 
the double wedge airfoil as a practical shape, but a bi- 
convex lens shape may be better structurally and at 
subsonic speeds. A rounded leading edge is definitely 
ruled out. The aerodynamic center is at mid-chord, a 
serious shift from the subsonic quarter-chord position. 

In comparing this theory with a more exact approach, 
the disturbance flow on either side of the wing may be 
regarded as a progressive wave traveling backward and 
outward from the leading edge. As in the unsteady 
one-dimensional case there is a simple general solution 
for this problem if entropy variations are ignored. As 
before, there is a functional relation between the veloc- 
ity variables; thus, if @ is the angle made by the flow 
with the x-axis and yu is the local Mach angle, then in a 
progressive wave,® * % 





6 = \ tan~'(A tan uw) — w + constant 
M2 = (y — 1)/(y + 1) (7) 
» = fly — x tan (u + 8)) 


The velocities are constant on a system of straight lines 
called Mach lines, each at the Mach angle to the stream- 
lines. For airfoils 0 is known on the surface, so, if the 


II 
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constant in Eq. (7) can be determined, yu is also known 
and the Mach lines can be drawn making this angle 
with the surface at each point (Fig. 1). At a discon- 
tinuity of the surface slope a number of these lines 
will be concurrent, since a range of values of « must be 
covered at this one point where @ changes discontinu- 
ously. With falling 0, u falls also, and the speed rises, 
as at the corner shown in the figure. At a corner con- 
cave to the flow, successive Mach lines would be further 
back so that each point would have two speeds, an im- 
possible state of things, and we can be sure that a shock 
wave must supersede the isentropic flow. Even at a 
continuous concave surface as illustrated, the Mach 
lines converge, and beyond where they meet isentropic 
flow is impossible; shock waves produced there, how- 
ever, may not affect the surface pressures. 

It remains to consider the effect of the bow shock 
wave, for any rear shock waves cannot affect the surface 
pressures. Rankine’s equations yield that, if the airfoil 
is pointed, then from a certain Mach Number onward 
the shock wave will be attached to the leading edge and 
oblique to the stream.” © Its behavior depends on the 
angle through which the stream must be turned at the 
surface. The smaller the angle, the sooner it happens, 
which is desirable since behind a detached shock wave 
there is subsonic flow and a stagnation point with a 
much larger drag. If the angle is 4°, the shock wave 
becomes attached at 1/ = 1.22, when the pressure in- 
creases by a factor y = 1.39 at the shock wave. The 
entropy change is 


¢,[log y + y log (y + A*) — y log (1 + A*y)] 


which is 0.0015c, in this case. Thus, even for a large 
pressure change, the entropy change is remarkably 
small, and little can be lost in assuming isentropic flow. 
So we may take @ as zero when yu has its main stream 
value, and the whole solution is then determinate. It 
can be shown that the entropy change is of the order of 
the cube of the angles turned through and that this is 
the order of error in the above approximation. Neg- 
lecting terms of this order throughout, we obtain a 
quadratic instead of a linear approximation, which is 
due to Busemann,*! who, with Ackeret,'* pioneered this 
theory; this modifies the pressure distribution given by 
linearized theory but leaves the total lift and drag un- 


altered. The aerodynamic center moves forward, how- 
ever, to (say) 45 per cent chord. 

For airfoils like the double wedge, a completely exact 
theory is possible as Epstein showed,” since part of the 
bow shock wave is straight and the flow at the surface is 
isentropic; the constant in Eq. (7) must be modified, 
however, to make it exact at the leading edge. 

In no case are the shock waves from the leading and 
trailing edges parallel; they point somewhat away from 
one another. The Mach lines from the surface inter- 
sect them, however, and each time that this happens 
they curve in toward one another and become weaker. 
The manner in which they do so is easily found if 
entropy changes are ignored.® At large distances both 
shocks form approximately part of a single parabola 
with the airfoil as focus. Both are ultimately at the 
Mach angle to the stream and of zero strength. At 
distance R they are apart by a distance of order »/ R and 
their strength is of order 1/+/R. Hence, the total 
momentum transport away from the body is the same 
at alllevels. Far away more and more air is influenced 
less and less, but the total influence hardly changes, and 
the whole region lies in the Mach direction to the main 
stream. Lastly, I must observe that in the paper I am 
quoting from I stated that the downwash in the wake 
itself is of order the cube of the maximum angle of de- 
flection. But Kahane and Lees,** correcting me, have 
shown that it is smaller still and that fourth power 
should be read for cube. 

Busemann remarked* that by use of a particular bi- 
plane wing shock wave drag could be eliminated entirely 
at zero incidence and one particular supersonic Mach 
Number and would be small for other Mach Numbers 
close to this. Unfortunately, lift will be hardly more 
than for either wing used separately. The linear theory 
is easily applied here; the previous theory is unaltered 
except at points behind the bow wave of both wings; 
here the potential is the sum of the two progressive 
waves, and two functions f and g must be determined, 
but the boundary conditions are adequate. 


CHANNELS 


Again, the theory may be applied to supersonic flow 
in channels if the velocity varies little therein. In the 
case of a straight channel isentropic deviations from 
uniform flow are found to be possible which are not 
damped out downstream but which oscillate, for the 
potential can be expressed as a Fourier series of terms 
in cos (nry/h) if the channel is 0 <y < h, and to satisfy 
the linearized equation the coefficient of this can only be 
a trigonometrical function of mmz tan u/h. Ward® has 
obtained similar results for circular tunnels and com- 
pared them successfully with an experimental pressure 
distribution. These facts indicate that great care is 
needed to obtain uniform conditions in the working sec- 
tion of a supersonic tunnel. Again, a two-dimensional 
supersonic jet, emerging into still air at lower static 
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pressure, and, hence, at far lower total head so that the 
boundary is practically an isobar, is of periodic shape,”* 
expanding and then contracting as the expansion waves 
reflect from the free boundary as compression waves, 
then expanding again, etc. (see Fig. 2). In the linear 
theory the expansion arms are replaced by straight line 
discontinuities and the jet boundary is made of straight 
lines only. The jet of circular cross section can only be 
treated by linear theory; the boundary in this case was 
found by Ward (unpublished) and is not exactly peri- 
odie but (in the technical sense) almost periodic. 

When an airfoil is placed in a supersonic tunnel, the 
bow shock waves are reflected at the walls, but if the re- 
flected waves intersect behind the airfoil, there will be 
no tunnel interference. When a shock wave is too 
strong, however, or at too large an angle to the wall, 
normal reflection at any angle is impossible, and a com- 
plicated state of affairs occurs known as Mach reflection 
with a three-shock intersection in midstream. 


SLENDER Bopy OF REVOLUTION 


Next I wish to describe how certain difficulties in 
the application of the linearized theory to flow past 
bodies of revolution, which are continually being redis- 
covered, may be surmounted. If the equation of a body 
of revolution in cylindrical coordinates is r = 6R(z), 
where 6 is its fineness ratio, the boundary condition is 


0¢/dr = [U + (0¢/0z)]6R’(z) 
on r = 6R(z). An approximation to this, no more 
drastic than those involved in the derivation of the 
linearized equation, is 
0¢/Or = UbR'(z) 
on r = 0, On this assumption the problem is easily 
solved, and the value of ¢ near the body is about 


us R(z)R’(z) log r + f log a d{R(y)R'(y)]$ 
\ 0 2(2 — y) f 
(8) 


Thus, while 0¢/Or on the body is UéR’(z) and so of 
order 6, 0¢/0z is of order 6? or, rather, 6? log 6. But 
linearized theory neglects the squares of the disturbance 
velocities, in particular (0¢/0r)?, which is of the same 
order as 0¢/0z. Thus, the methods we have used are 
inconsistent.*! Again, in the pressure, given by 


dp | 1 o¢\? 3 (28) - 
ft +3(u+%) +3 = = constant (9) 


the third term capot be neglected when compared with 
Ud¢/dz. 

The correct answer to these anomalies is that the 
value of ¢ obtained by the method is right but that to 
obtain the pressure the approximation 


pb — po = —po[U(0d/d2) + (1/2)Q¢/dr)?] 


should be used. This statement is proved by using the 
exact equation of isentropic flow and expanding the 
potential in powers of 6. It is found that simple 
powers «re insufficient and that, as well as terms in 6°, 
54, 6°, ..., we must have terms in 6‘ log 4, 5° log 6, 5° log? 
6, etc. Each coefficient is a function of 7 and 2, and its 
behavior for small r can be found uniquely by use of the 
boundary condition and the condition that waves be 
propagated outward and backward at infinity. When 
this is done, a series for the pressure coefficient (p — 
po) /(1/2)poU? on the body is rigorously found, and the 
first terms were given by Broderick.*® They are of 
orders 6°, 5? log 6, 64, 64 log 6 and 64 log? 6, the largest 
term neglected being 0(6° log* 6). Now the terms in 6? 
and 6? log 6 are precisely those found by the corrected 
linearized procedure indicated above, which is therefore 
the right first approximation. The only gap in the 
argument is the assumption of isentropic flow, but the 
entropy change is of order the cube of the pressure 
change behind the shock and so presumably of order the 
cube of the pressure coefficient on the body—that is, 
0(6° log* 6), which is already neglected. But to put all 
beyond dispute the pressure change at the shock wave 
for cones of small fineness-ratio, 6 has been investigated 
analytically by two authors® * using the exact equa- 
tions with the change in entropy, and it is found to be 
(3/2) y(y + 1)M*®(M? — 1)—'po64 for small 6, varying as 
the fourth power of the fineness-ratio—an extremely 
weak wave—and the theory checks Broderick’s form of 
the pressure coefficient. This has also been checked 
against various exact solutions for the cone found 
numerically” at M.1.T., following Taylor and Maccoll ;* 
it is close indeed for cones up to 10° semiangle, while the 
correct linearized terms give a better approximation 
than do other approximations, such as the exact solution 
of the linearized equation with exact boundary condi- 
tions which von Karman,”* to whom most of the credit 
for this work is due, originally advocated. 

Similar work has been done for the body of revolution 
at yaw,*‘ when an additional expansion in powers of the 
angle of yaw a@ is required. Terms in the pressure 
coefficient of orders ad, aé*, aé® log 6, a®, and a*§ have 
been found in addition to those mentioned previously. 
The drag on a body of revolution, being produeed by 
pressures of order 6? on an area of order /*5*, where / is 
the length, must be proportional to pU*5‘/*; I will give 
its precise form later. If it has a flat base, there is, in 
addition, a base drag due to suction by the dead air 
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behind the base. At yaw there is an additional drag 
(1/2)pU*a? times the base area, and the lift is pUa 
times the base area as Tsien found*® in 1938. The terms 
in a*® in the lift are found to be insignificant, but the 
terms in a6’ are larger and cause a fall in lift with Mach 
Number for thin bodies and a rise for thicker ones. 
The aerodynamic center is at the center of mass of the 
air displaced by the body from a circular cylinder drawn 
through the base. For a body with a pointed tail, 
however, the lift is much smaller, and the most powerful 
force is a pure yawing couple. 


THREE-DIMENSIONAL SLENDER BODIES 


Ward® has generalized the linearized theory to cover 
the flow past any three-dimensional body slender in all 
directions perpendicular to the stream but not neces- 
sarily of revolution, I will give a short version of this, 
remarking that the bold steps therein are justified by 
the careful work previously described and by similar 
work, The equation 

Op  O'h _ 0p 


ae ap * ng 


we = p*cot? yd = (10) 
in Heaviside’s operational notation is solved in cylin- 
drical coordinates (r, 0, 2) as 


—Ko(pr cot u)fo(z) + Kilpr cot we“fi(z) + ... 


where the Bessel functions K,,** not J, must be used be- 
cause of the form of the waves at infinity. For small r 
this is approximately 


[log (1/2) pr cot u + C]fo(z) + (pr cot u)~'e"fi(z) + 
2(pr cot pw) —29218 8) cs 


This is a harmonic function of x and y expanded in 
powers of r, beginning fo(z) log r + k(z) + ..., where 


k(z) = [log (1/2)p cot » + C]fo(z) = 
- cot yu u , 
I Bs As — nh (t)dt (11) 


The boundary condition is 0¢/0n = Udn/dz, where, as 
Fig. 3 shows, dn is the distance along the normal be- 
tween sections of the body at a distance dz apart. 
This boundary condition determines a harmonic func- 
tion to within a constant term, but luckily the constant 
term is already known by Eq. (11). Further, the 
coefficient of log r, or “‘source strength,’’ is 





U dS (12 
a= zZ 
Qr dz 2Qxr dz 


1 edo,  U yedn 


s 
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where d.S/dz is the rate of change of cross-sectional area 
in the z-direction. The drag, under certain restric- 
tions, is obtained quickly in terms of the momentum 
transfer across a small cylinder as 


l dg OS f 
Sen 2 — a —2? 0 Z) 
f wal Oz . 7 0 pu & 


[fo’(s) logr + k’(z)|dz = 2mp f fo'(s)k(s)ds = 
0 


er 1 
. fi s*(eds f log Raia as? S” (t)dt (13) 
27 J 0 0 s—i 


which extends von Karman's form for bodies of revolu- 
tion to slender pointed bodies of any shape. Other 
total forces are obtained by Ward® under various condi- 
tions, including the presence of wings of small aspect 
ratio for which the above theory must be modified. 


QUASICYLINDRICAL Ducts 


In 1945, with an eye on the athodyd engine, I gave®! 
a theory of supersonic flow past ducts of a shape which I 
called quasicylindrical—that is, with radius not varying 
greatly from a fixed value R. The form of the pressure 
coefficient is 


, zs—s \ n(s)ds ; 
t 3) — f w( —) as] (14) 
-" on 0 R cot u/ R cot u ii 


where n(z) is the angle made by the outer surface with 
the main stream and WV is a new function that was tab- 
ulated. The n(z) term corresponds to the two-dimen- 
sional airfoil solution, but the integral grows as we re- 
cede from the front and finally dominates the solution. 
Later,® I used this theory and the properties of the 
function W to get over the deficiency of the theory of 
thin bodies of revolution which I have been discussing: 
when no discontinuities of shape were allowed on the 
surface. It breaks down when these are present (as for 
a double cone) and gives an infinite drag. The correct 
theory is found to give a drag rather larger than when 
the surface is smooth, and the drag coefficient depends 
on Mach Number, which it did not before. The large 
pressure drop behind a discontinuity in slope is made up 
rapidly, which may adversely affect the boundary layer. 
This theory also applies to ducts of any shape if they 
are thin compared with their total length. 

Ward* improved the theory of quasicylinders and ex- 
tended it to cover the flow inside, for which my theory 
was incorrect; but here again linearized theory steers 
through difficult waters. The disturbance propagated 
from the lip down the Mach cone, which is a discon- 
tinuity in velocity or velocity gradient, is reflected at 
the axis into an opposite-pointing cone on which the 
disturbance is not a discontinuity but a logarithmic in- 
finity. On reflection at the surface and second reflec- 
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tion at the axis, it becomes a discontinuity again but of 
opposite sign. The problem of using the nonlinear 
equations to find what really happens on these reflected 
characteristics is one of the most interesting unsolved 
problems in high-speed gas flow but is even harder than 
at first appears. It can only be said that the logarith- 
mic infinity is actually impossible on a characteristic 
and that the series of which the linearized solution is the 
first term diverges near it; but the solution may be a 
reasonable approximation elsewhere. The problem is 
important, since the singularity presents a barrier to the 
use of the method of characteristics. 


LINEARIZED THEORY OF THREE-DIMENSIONAL WINGS 


Next I must speak of the supersonic flow past three- 
dimensional wings and offer some guidance, however 
presumptuous, through the mountains of literature on 
this topic. Treated as a linear problem it can be split 
up into one symmetrical about the plane of the wing 
and another antisymmetrical—that is, into the problem 
of the lift on a flat wing at incidence with given plan 
form. Some geometrical considerations are often help- 
ful. Thus for a wing of rectangular plan form the flow 
over parts of the wing will be identical to the two- 
dimensional airfoil flow, for a point is only affected by 
conditions in its forward Mach cone and these are the 
same in the two flows. Thus only the parts of the 
wing in the Mach cones from the tips will feel end 
effect. 

The drag problem is mathematically simple. With 
the main stream in the z-direction and the wing in the 
(x, ) plane with equation y = Y(x, 2) we have approx- 
imately 0¢/0y = UOY/dz on the wing surface, and so 
also to a first approximation on the central plane of the 
wing, if @ is the potential above the surface. On the 
rest of the plane y = 0, 0¢/Oy vanishes, and so we have 
an equation (0¢/Oy),-0 = f(x, 2) everywhere. To 
solve the linearized equation with this boundary condi- 
tion is straightforward. Thus a solution of Eq. (10) is 


f A(n,p)e™e to" cot MM dy 


fo) 


where the negative exponential must be taken to ensure 
that wave propagation at infinity is in the right direc- 
tion. By satisfying the boundary conditions and inter- 
preting the operational form we deduce”* that 


1 *x+2 tan yw s—|x—é|cot wp 
fa fw x 
TJ x—s tan “ 0 


(¢),=0 ee 
du 


Vi(z — u)? — (x — 2)? cot? yp] 





(15) 





The drag is then 


wf (5), clos) 
-2 ix ~ =) @ 16 
rf % ¥ 0 (* y=0 Oz y=0 “ ( ) 


to evaluate which only patience is necessary. 


The drag on a rectangular wing of aspect ratio greater 
than tan y is found to be the same as the value given by 
two-dimensional theory, For lower aspect ratios the 
Mach cone from each tip intersects the other side of the 
wing and the drag is lower.”* But we may adopt a 
policy of placing the leading edge at an angle not (1/2)x 
but (1/2)" — 7 tothe stream. 7 is then called the angle 
of sweepback. The advantage of this for an infinite 
wing follows from the principle of relativity, since, to an 
observer traveling along the wing with velocity U sin r, 
the flow will appear to be a flow with velocity U cos r 
normal to the leading edge, so that transonic troubles 
will arise only when U cos r approaches the velocity of 
sound. When it exceeds it—i.e., JJ > sec r—the drag 
per unit leagth is 


cos? r VM? — 1/V iar — sec? 7 


This ex- 
Practical 


times the drag of the same wing not vawed. 
ceeds unity only in a small range of 7. 
shapes embodying this idea of sweepback are the arrow- 
head wing and the sturdy delta wing (Fig.4). For these 
the drag, even on this theory, is finite at 1/ = sec r. 

It is simplest for comparing different plan forms to 
assume that each section parallel to the stream is, say, 
a double wedge. In fact, wherever this section is a 
polygon, the expression for the drag is simplified, since 
f(x,z) is a step function of z, and in differentiating 
(¢),=9 with respect to z, only terms due to the dis- 
continuities of f contribute to the result. If these are 
atz = a(x), 2 = (x), etc., and are of magnitude 
a(x), ao(x), etc., we have 


Pa) 1 x+ztanyp 
Bh dn * 
Oz y=0 T x—stany 


a, off) dt (17) 
[z — c,(t)]? — (x — £)* cot? u} 


the sum being for values of r with z — c,(t) > |x - t| cot p. 
The drag becomes 


-¥ p Be & J / an(x)a,(t) X 


cosh~! eee dxdt (18) 
is = t| cot u 








summed and integrated over values of n, r, x, and ¢ such 
that the quantity in brackets exceeds 1. One imme- 
diate conclusion is that drag is unchanged if the flow past 
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a wing is reversed, a fact that often reduces labor. For 
a double wedge section there are only three values for n 


or for r. Thus for a delta wing of maximum chord c, 
we have 
! 
ex(x) = |x| tan 7 
/ | 
c(x) = (1/2)(|x| tan r + c) 
s(x) = 


while if the thickness-ratio of the sections is 6, 


Gx) = U6, aa(x) = —2U6, Galt) = U8. 


The ratio of the drag to the drag of a rectangular wing 
with the same area and shape of section,® which would 
be 2p U6" tan uc? cot 7, depends only on 


n = cot rcot wu = cotr VM? — | 


Its behavior is shown in Fig. 5; it never exceeds 1.24 
and is greatest when the Mach cone lies along the lead- 
ing edge. As » — 0—i.e., J — 1—the function be- 
haves like (2/7)n log (e/n) and so the drag is like 

: pU5°c? cot? r log nee = 

7 V (M? — 1) 
which becomes large as JJ — 1. This is a much milder 
singularity than that occurring for a two-dimensional 
wing and is of the same order as that for a body of 
revolution whose surface has discontinuities in slope. 
We must expect that the theory will be incorrect near 
M=1. , 

Squire® suggested a system of wing sections for a 
delta wing with rounded leading edges if the wing is to 
be flown at a Mach Number less than sec 7, since the 
flow normal to the leading edge will then be of subsonic 
type and no trouble such as a detached shock wave 


should result. With a wing surface given by 


y = (1 — |z/cl) W[(¢/c)? — (x tan 7/c)?] 


- 


the central section is a biconvex arc of thickness /, and 
all other sections are of the shape of a conventional air- 
foil of ‘‘Karman-Trefftz”’ type. The ratio of the drag 
to that of a rectangular wing of the same area with 
section similar to the central section, which would be 
(1/2)pU*(16/3)¢? cot 7 tan yw, is shown in Fig. 6 as a func- 
tion of 7 in the range n < 1—i.e., MJ < sec r—in which 


the wing can be used. The rounded leading edge will 
at least be an advantage at subsonic speeds. 


BUSEMANN’S CONE-FIELD THEORY 


In the lift problem with the flow antisymmetrical, we 
have ¢ = 0 when y = O where the wing is not, but 
(0¢/Oy), -9 is given as before on the wing surface. This 
mixed boundary value problem is much harder, and no 
complete solution in a closed form has been given. The 
literature consists of a number of special methods, each 
suitable for particular classes of problems. The most 
developed is Busemann’s cone-field theory.®! This is 
applicable to problems where the velocity vector is con- 
stant on every straight line through some particular 
point. This is true in the lift problem for a delta wing, 
the point being the apex, since neither equation nor 
boundary conditions are altered if a change of scale is 
made. It is true also inside the Mach cones from the 
tips of, say, a rectangular wing and even in the region 
where these intersect the flow can be obtained by adding 
up these two flows and subtracting the two-dimensional 
flow. 

In a cone-field problem the velocities are functions of 
x/ztan uw = x, and y/stanyu = y, alone. In the (x, y;) 
plane the Mach cone is the circle x; + y,> =1. Ifx, = 
r cos 0, y = rsin 6, and we put p = (1 — V1 — 7r’)/r, 
which increases from 0 to 1 as r does, then each velocity 
component, say #, is found to satisfy Laplace’s equation 


(07u/Op*) + p~'(Ou/Op) + p—*(02u/00?) = 0 


so that if ¢ = pe”, u, v, w are the real parts of regular 
functions of ¢ say U(g), V(g), W(e). These are not 
independent but satisfy U’(f): V’(O): W'(O) = —-(e+ 
1): 4(¢? — 1): 2¢tany; so that there is a sort of general- 
ized potential F(¢) such that 


Us) = — S (8+ 1) Fd 
Vis) =i S (8 — 1) Fd (19) 
Wr) = 2tanu Sf ¢F(e)de 


To determine /(¢) from the boundary conditions is now 
a matter for conformal transformation theory. The 
conditions apply on the Mach cone and at the wing— 
that is, on the unit circle and the real axis. Ward and 
Goldstein® show how normally a transformation of 
the upper semicircle into the whole upper half plane 
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gives F(¢) with a minimum of effort, the integration to 
obtain u, v, w being left to the end. This integration 
always produces elliptic integrals, but, if it is left to the 
last, the main body of the work is performed by ele- 
mentary reasoning that gives the method an advantage 
over some versions, in which elliptic integrals and their 
transformation theory predominate from the start. 

The cone-field theory can be extended to cases where 
the velocities are not independent of radius but vary as 
some power of it, which is specially useful in calculating 
aerodynamic derivatives. Other problems can be 
tackled by superposing a continuous distribution of 
cone fields’* or by a general method for large enough 
aspect ratio due to Evvard,”! which uses Eq. (15). As 
the aspect ratio becomes small compared with tan yg, 
the problem becomes progressively more complicated, 
as is familiar from transient diffraction problems in the 
theory of sound with some of which the present theory 
is identical, an analogy that Gunn’? used to solve prob- 
lems where the aspect ratio was less than tan yp. 

Passing to results, I should say, first, that for a rec- 
tangular wing the lift falls to zero in the parts of the 
wing inside the Mach zone from the tips in such a way 
that the total in these parts is half the two-dimensional 
lift. These two triangles each have area (1/2)c? tan p, 
so with span > the lift coefficient is 

4a tan p[1] (1/2)c? tan (u/bc)]} 
This is correct only if ¢ tan uw» < 6. When b/c has 
dropped to (1/2) tan yw, the factor has dropped to 1/5 
Results are also known for trapezoidal shapes; if the tip 
triangles are cut out altogether, the lift coefficient takes 
its two-dimensional value, but this plan of Busemann’s 
is probably inferior on the whole to his other idea of the 
sweptback wing. The lift coefficient on a delta wing® 
is 4a tan uw when JJ > sec r—i.e., 9 > 1—and is 27a 
cot r/E’(n) when n < 1, where £’ is one of the elliptic in- 
tegrals satisfying £’(0) = 1, E’(1) = (1/2). In cases 
where any part of the trailing edge is at an angle greater 
than the Mach angle to the main stream, the solution to 
the lift problem is not unique unless Joukowski’s condi- 
tion of finite velocity be applied over this portion of the 
trailing edge. The drag due to incidence equals the lift 
multiplied by the angle of incidence (since all pressures 
are normal to the plate), unless part of the leading edge 
is at an angle greater than the Mach angle, when the 
infinite suction there, familiar and known to give correct 
results in subsonic theory, will reduce the drag due to 


incidence. 


UNSTEADY PROBLEMS 


A number of unsteady problems, such as the effect on 
a wing’s motion of small amounts of pitch, roll, yaw or 
sideslip, are reducible to steady problems by techniques 
familiar in low-speed flow—that is, by their reduction to 
the calculation of aerodynamic derivatives, which are 
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found by varying the angle of incidence of different 
parts of the wing. The genuine linearized equation of 
unsteady flow is 

0*¢ yo? 0*6 


We epi 9 
“oa a CO 





studied in the two-dimensional case by Possio,*® von 
Borbely,* and by Temple and Jahn.*° Using solutions 
built up from a fundamental one 


eft U “is sect) F51 lw sec? uw X 


V (5? sin? sin? , Me 1 y? y? cos? s? n)] (21) 


they describe fully the behavior of a harmonically 
oscillating airfoil at supersonic speeds, which is im- 
portant because of the danger of aerodynamic flutter 
and which is extended to the three-dimensional case by 
Garrick and Rubinow.” 


NONLINEAR TRANSONIC APPROXIMATION 


An approximation to steady transonic flow that is 
partly subsonic and partly supersonic, differing from 
the linearized one in one point only which, however, 
makes -it nonlinear, has been much discussed re- 
cently; 88 but its value is uncertain, since velocity 
variations are known to be great (however thin the 
body) at transonic speeds. But the inference it affords, 
that the critical Mach Number at which a precipitate 
drag rise begins is less than 1 by a quantity roughly pro- 
portional to the two-thirds power of the thickness, is 


probably correct. 


THEORY OF CHARACTERISTICS 


The theory of characteristics provides a convenient 
numerical method of calculating particular gas flows. 
It is applicable only when the flow has two independent 
variables, such as (x, y), (7, 2), (x, ¢) or (7, t), and when 
disturbances are propagated at a finite rate, which ex- 
cludes steady subsonic flow in which a change of bound- 
ary affects all parts of the field at once. 

For a partial differential equation of the first order, 
all solutions and all their derivatives are continuous 
everywhere where this is true of the coefficients, but it is 
not so with second-order equations. Where these are 
linear in the second-order derivatives,*‘ as in 

076 p22? 
Ox? “| Oy 
where A, B, C, D depend on x, y, ¢, 0¢/0x, and 0¢/Oy, 
solutions are possible with discontinuities in some of 


their derivatives across certain lines. If dx,dy is a 
small displacement along one such line, then 


Ady? — Bdxdy + Cdx*? = 0 


For a given solution, lines satisfying this condition are 
called characteristics; they are a fixed set of curves only 
if A, B, C depend only on x, y._ If four arcs of charac- 


+D=0 (22) 
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teristics form a quadrilateral PQRS, the solution within 
is determinate when known on PQ and PS. If PQRS 
is small and the solution is known at P, Q, and S, we 
may interpolate for its values in between and obtain the 
value at R and also the position of R. Now, if given 
any three points the quadrilateral may be completed, it 
is clearly possible to use this to obtain progressively the 
solution in a whole field, starting from the solution at a 
number of points at one end; boundary conditions may 
be satisfied on the way. This delightful progressive 
method is in contrast with the “‘relaxation’’ methods 
necessary in steady subsonic flow (where the character- 
istics are unfortunately imaginary) in which continued 
adjustments over the whole field are necessary in a 
desperate effort to satisfy all the boundary conditions 
simultaneously. 

In plane isentropic supersonic flow, the equation of 
motion is linear if «, v (the velocity components) are 
taken as independent variables, and the characteristics 
are fixed curves in the (7, v) plane; Prandtl and Buse- 
mann” used this to expound a rapid reckoning process 
useful in designing supersonic tunnels. In axisym- 
metrical flow the labor is greater, and, when entropy 
variations are allowed, the situation is much worse be- 
cause discontinuous entropy gradients can be propa- 
gated along streamlines which are therefore a third fam- 
ily of characteristics. But Crocco,** Tollmien,** Gud- 
erley*' 74 (especially), and (more recently) Meyer* 
have done their best for us in this case. 


HopoGRAPH METHOD 


I have left until the last the hodograph method, in 
which the linearity of the equations of plane steady 
isentropic flow, subsonic or supersonic, when 1%, v, or 
better g, 0 (the magnitude and direction of velocity) are 
taken as independent variables, is systematically used. 
This method has been regarded as man’s white hope of 
understanding transonic problems other than by flight 
tests or by comprehensive application of numerical 
procedures based on a combination of relaxation and 
characteristic techniques. However, it is applicable 
only to steady two-dimensional flow where vorticity 
and entropy variations can be neglected, and, further, 
it seems itself to encounter serious difficulties in the 
transonic region. But these difficulties may be sur- 
mountable, while those of the approximate methods are 
inherent. 

Using a stream function y, those value at each point 
measures the rate of mass flow between the streamline 
through that point and some fixed streamline, so that 





oy /Oy = pt, oy /Ox = —pv 
we can easily obtain the equation for y in terms of ¢,4 as 
0 /q ~) o*y 
—{=-—]) = (4-1 2 
pq 5 (2 i ( ) v7 (23) 


That for ¢ is similar though slightly less convenient. 


Writing Bernoulli’s equation in the form 
(1/2)q? + [a*/(y — 1)] = (1/2)q?mar. 


so that Qmar, is the maximum velocity attainable in the 





field, Chaplygin® gave as a solution y = y,(q)e’”, where 
¥n(Q) = Q"F (Anda; m + 1; 9?/Q*maz.) 
1 n(n + 1){ (24) 
a, +b, =n — , a), = — ——— 
<9 2(y — 1) 


Here F is a hypergeometric series, being equal to 1 when 
q/Qmar, = O—i.e., at zero Mach Number—when the 
solution becomes the familiar solution g"e'’. Chap- 
lygin showed also that, if a solution of Eq. (23) is 
obtained, to obtain the points in the physical plane 
corresponding to different values of g, (as is necessary 
for interpreting the solution) the formula 


° Sie ae 1 oY  *) 10 95) 
x+tw= ie (x +5 50 e”dé (25) 


may be used. 

Chaplygin extended to gas flow by these means the 
solutions of free streamline problems given for low-speed 
flow by Helmholtz, Kirchhoff,* and others.*” Thus, in 
the problem of the efflux of gas from a slit in a plane 
wall, we assume constant pressure and, hence, constant 
speed on the jet boundary so that the boundaries of the 
flow are gq = U (say) and @ = +(1/2)z. For incom- 
pressible flow the problem is easily solved by con- 
formal transformation in terms of the complex variable 
t = ge~"/Uas 


¥ = k&[log ¢ — log(1 — ¢?)] = 
—k (0 oe fou 2nd) (26) 


i na U™* 


The constant k can be found in terms of the width of the 
slit. Chaplygin showed that in gas flow the correspond- 
ing solution is 


1 Yon(q) ‘ ) al 
= —k \0 — ———— sin Sad 27) 
y ( -+- Ms te, (U) sin 2n (27) 


It is in fact clear that this will take a constant value on 
6 = +(1/2)z and that also, when g = U, Eqs. (26) and 
(27) are identical. Again k is found in terms of the 
width of the slit using Eq. (25). The ratio of width of 
jet to width of slit deduced from this is shown, as a 
function of the Mach Number on the jet boundary, in 
Fig. 7. The method applies only when this Mach 
Number is less than 1; otherwise Eq. (27) does not con- 
verge. In this way any free streamline problem soluble 
for low-speed flow can be extended to subsonic high- 
speed flow with ease if a series expansion like Eq. (26) 
can be found valid in the whole field, and by special 
methods in some other cases, but no solution has been 
found with even a local supersonic region. Chaplygin 
also gave an approximate form to the hodograph equa- 
tions, and other authors have since given better 
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ones,” 4° 44 49 which are handy for getting rough esti- 
mates of the effects of compressibility. However, I 
shall say nothing of these because they constitute no 
theoretical advance on the approximate methods I 
have mentioned earlier. 

It was soon realized that, by finite combinations of 
solutions y,(g)e’”, an indefinite number of exact steady 
irrotational gas flows could be constructed with ease, 
which appeared a tremendous advance. However, it 
later became clear that they were not valuable, since the 
solutions, though without singularity in the hodograph 
plane, had several in the physical plane; further, they 
all had zero velocity at infinity. The most instructive 
was that due to Temple and Yarwood™ (following 
von Ringleb**) who took 


¥ = Ag™[l — (g°/g°maz.)]"'"~” sin 0 


In the physical plane the flow can be regarded as flow 
through a channel, in which the speed rises from zero at 
one end to a maximum value and falls to zero at the 
other end, the stream moving through two right angles 
in the process. Iu Fig. 8 any streamline further from 
AB than CD could be chosen as boundary with no 
qualitative effect. But none beyond AB can be chosen, 
for each such streamline if calculated behaves madly 


when it reaches each of the cusped lines shown. It 
passes over the first branch, turns back at the second 
and turns forward again at the first. The flow above 
AB therefore is physically impossible, and the cusped 
line whose appearance has limited the field of flow is 
called a limit line. It is a mathematical obstacle to 
isentropic flow taking place under certain conditions; 
now a shock wave is a physical obstacle to the same 
thing. This correspondence had been clear for some 
time when Guderley”* showed that, given a flow (so- 
called) with a limit line, it is mathematically possible to 
find a unique position of a shock wave starting at the 
cusp of the limit line and continuing not along either 
branch but usually in between, so that the anisentropic 
flow behind the shock wave is continuous with the isen- 
tropic flow below the streamline through the cusp. 
The practical calculations are prohibitively long but the 
fact is illuminating. 

The hodograph method was applied in a more 
advanced manner to the flow in supersonic wind tunnels 
by Frankl® and myself” independently. A new sing- 
ular curve was found called a “‘branch line,’’ which is 
almost the opposite of a limit line since in this case the 
physical plane is nonsingular but the hodograph plane 
is singular. The branch line has a cusp in the hodo- 
graph plane, but not in the physical plane where it 
touches the sonic line once but lies in the supersonic 
region. We invented separate methods of getting 
round the singularity which makes the work a possible 
basis of a method for designing tunnel shapes. 

Craggs showed*! that the only possible singularities 
in the hodograph transformation are (a) isolated ones 
in the subsonic region, (b) limit lines and branch lines 
in the supersonic region. A limit line is an envelope of 
isobars and of one family of characteristics, while the 
streamlines and the other characteristics have cusps on 
them, except at their own cusps through which the latter 
curves proceed regularly but with infinite curvature. 
A branch line is itself a characteristic on which isobars 
touch characteristics of the other family, but it does not 
envelop these. The condition ‘‘isobars touch charac- 
teristics’ is thus not a test for either. 

Finally, the method was applied to flow of a uniform 
stream round a body. At first it was thought that 
Chaplygin’s method could be used—that is: express ¥ 
for incompressible flow as a series of terms (¢/U)" sin n@ 
and replace these by 


[¥n(Q)/¥n(U)] sin 16 (28) 


throughout. But there were two objections: 

(a) Since the boundary no longer consists of curves 
q = U and @ = constant, there is no reason why this 
solution should satisfy the boundary conditions. 

(b) Since y is infinite when g = U and 6 = 0, it is 
impossible that one such series can represent y every- 
where; thus we will have an ascending power series 
when g < U and a descending one when g > U, and, 
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since the two parts of the plane with g < U are discon- 
nected, at least four series will always be necessary to 
describe the whole flow, each valid in a different region. 
Now, there is no reason why the modified series taken 
together will represent a continuous flow. 

The result of (a) was, firstly, that the use of the cum- 
bersome denominator y,(U) was pointless and any fac- 
tor could be inserted which would make the series con- 
verge in the appropriate region, and, secondly, that the 
hope of finding the flow past a given body was aban- 
doned in favor of that of finding some exact flow past 
any body or, less vaguely, a flow reducing for zero Mach 
Number to a flow past a given body. The result of (b) 
was that more elaborate methods had to be used which 
would give a single expression for y valid throughout 
the plane, from which series in y,(qg) sin 7@ might be 
obtained in different parts of the plane if they proved 
easier to compute as they often do. But these series 
are not simply the incompressible series modified as in 
Eq. (28) but contain additional terms. Several in- 
correct theories were produced, but those of Berg- 
man,’ Cherry,” and myself are correct. The last 
two are effectively the same, and they show that 
special methods have to be used when there is circula- 
tion round the body, when the factor 1/y,(U) in Eq. 
(28) must be replaced (in a certain sense) by 


fy_.(U) + Upl,(U)]/ — n) 
In the absence of circulation a much simpler factor 
than either is permissible. 

When a solution is computed by one of these methods 
for different values of LU’, the body shape varies some- 
what with Mach Number. When J/ reaches a ‘“‘lower 
critical’ value,” a sonic speed first appears in the flow; 
when J/ exceeds a “higher critical’? value (still sub- 
sonic), there is a limit line in the flow which has seemed 
to pierce the surface from the inside, cusp first. Buta 
warning is necessary: ‘The body shape changes with J/ 
and has an infinite curvature at the higher critical value, 
so no meaning can be attached to the “higher critical 
Mach Number for a given body”’ with finite curvature. 
Probably for such a body, as the Mach Number in- 
creases, the point at which isentropic flow becomes im- 
possible and a shock wave becomes necessary is not 
characterized by any singularities in the flow field. 
Thus my own personal view is that limit lines will not 
help us to understand the flow about airfoils but that 
some other development of the corrected hodograph 


theory is required. 
CONCLUSION 


In conclusion, I will say that our understanding of the 
high-speed flow of gases is growing rapidly; as a mathe- 
matician I believe that our endeavors in this field are no 
waste of effort or of mathematical technique not only 
because we are assisting in a great new engineering 
‘but also since we are 
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getting to grips in these problems with that old bogey. 
man, the nonlinear partial differential equation, and 
smelling out his ways in a manner for which our col. 
leagues, in the more fundamental parts of physics, may 


later be grateful. 
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On the Addition of Heat to a Gas Flowing in 
a Pipe at Subsonic Speed’ 


JOSEPH V. FOAT ann GEORGE RUDINGER?t 


Cornell Aeronautical Laboratory 


ABSTRACT 


It is shown that, contrary to widely accepted opinions: 

(a) The addition of heat to a subsonic flow in a pipe modifies 
the flow conditions both downstream and upstream of the region 
of heat addition. 

(b) The Mach Number of the heated flow does not necessarily 
increase, approaching the sonic state, when heat is added; _in- 
stead, it will increase or decrease, depending on the boundary 
conditions. 

(c) The static temperature of the gas does not decrease when 
heat is added within the range of initial Mach Numbers from 
1/+/y to 1.0; instead, it increases continuously. 

(d) The attainment of sonic velocity in the heated region does 
not limit the amount of heat that can be added. 

(e) The mass flow may decrease as a result of the addition of 
heat but this decrease is in no way related to the attainment of 
sonic velocity in the heated region. 

A quantitative treatment of the problem is presented for the 
cases of greatest practical interest, and numerical examples are 
given to illustrate the procedures and conclusions. 


INTRODUCTION 


T IS OFTEN OF INTEREST to determine how a given flow 
of a gas in a pipe is modified by the addition of heat. 
This problem has received increasing attention in re- 
cent years. 

Villey’ derives the logarithmic derivatives for the 
state parameters and velocity in terms of Mach Num- 
ber and heat added and concludes that in the range of 
Mach Numbers between 1/+/y and 1.0 there occurs an 
inversion of the thermal effect of heat addition or, in 
other words, “‘a process with negative specific heat.”’ 

Similar conclusions are reached by Focke,? who as- 
sumes that the state at the inlet be the known boundary 
condition. As a result of this he concludes that the 
condition .\J = 1 sets a definite limit on the amount of 
heat that can be added. 

Hicks,* and Hicks, Montgomery, and Wasserman’? 
treat the problem in differential form and reach the fol- 
lowing conclusions: (a) As heat is added to a subsonic, 
compressible flow, the Mach Number will always in- 
crease, and (b) the temperature increases only when 
heat is added at M < 1/+/y but falls in the range of 
Mach Numbers between 1/V¥ and 1.0. These changes 
are properly computed and interpreted as gradients 
within the heating region in a steady régime, and the 
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flow modifications that result from changes of heat 
addition are not considered. 

Szezeniowski‘ again derives and integrates the dif- 
ferential equations describing this phenomenon, assum- 
ing a fixed inlet state as the known boundary condition. 
He concludes that “‘if the gas has an initial velocity 
equal to the velocity of sound, then heating is out of the 
question”’ and even suspects that a-gas could not be 
heated or cooled between JJ = 1/Vy and M = 1.0— 
i.e., that the heat conductance coefficient would equal 
zero between these limits. He expresses surprise at 
these theoretical results and recommends experimental 
verification of his conclusions. 

Hawthorne® '® and Shapiro'® find that, for given 
initial values of Mach Number and temperature, there 
is a maximum possible rise in stagnation temperature— 
i.e., a limit to the amount of heat that can beadded—and 
suspect that the flow would be choked if the stagnation 
temperature were increased above this limiting value. 

Chambré’ and Lin” ° also study the problem in dif- 
ferential form and conclude that, ‘‘when heat is being 
continually added, a subsonic flow cannot become super- 
sonic but will always approach the sonic state.’’ The 
singularity of the logarithmic derivatives of the state 
parameters at 1J = 1.0 means, in their opinion, that 
“either the gas is unable to absorb heat, or rather the 
assumed conditions cease to exist (e.g., combustion 
stops).’”’ The mechanism by which combustion would 
stop is explained in terms of a “‘sharp flame front,” 
which would be established as soon as the Mach Num- 
ber reaches 1.0. They also state that ‘‘the temperature 
actually decreases when heat is being added if 1/-/ 7 < 
M <1.” This conclusion is also reached by Wood.’ 

Bailey" derives a relation for the pressure gradient 
which may, in his opinion, present singularities that 
would explain the often observed “rough burning.” 
Fulton!’ challenges this view and finds that the singular- 
ity of the pressure gradient only occurs at MJ = 1.0, 
whereas many combustors in which J/ never reached 
1.0 have exhibited rough burning. Fulton also states 
that the addition of heat at constant area between M = 
1/o/y and M = 1.0 results in a decrease of static tem- 
perature. 

Roy'® arrives independently at a conclusion that is 
similar to those of Szczeniowski and Chambré and Lin: 
at 1 = 1 no heat can be exchanged with the flow and 
no chemical reaction can take place. He also inter- 
prets his finding as indicating that the combustion that 
may take place in a flowing mixture becomes weaker 
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and weaker as the initial Mach Number approaches 
1.0, until the flame may be blown away or assume a 
pulsating character. He finds that practically no com- 
bustion could take place in a mixture that is ignited at 
M > 0.6 and recommends that the Mach Number be 
kept as low as possible in the region of ignition. Gener- 
ally similar conclusions are reached by Samaras; 
Graham;'* Turner, Addie, and Zimmerman;!” Zu- 
crow;'® and Barry.” 

A different view is taken by Kahane and Lees,'® who 
demonstrate the fallacy of the assumption of fixed inlet 
conditions, through an analysis of the transient that 
follows a sudden addition of heat. This analysis is 
carried out for two specific cases, and it is shown that 
the transient leads to considerable modifications of the 
inlet conditions in both cases. 

The problem of heat addition to a gas flowing at 
subsonic speed in a duct of constant cross-sectional 
area is re-examined in this paper. It is shown that 
the addition of heat modifies the flow parameters not 
only downstream but also upstream of the region of 
heat addition, in a manner that depends on the bound- 
ary conditions. As a result of this, the conditions up- 
stream of the region of heat addition should not be 
taken as the initial conditions (before heating), and the 
equations derived in previous publications can only be 
interpreted as relating the conditions upstream and 
downstream of the heating region and not before and 
after the addition of heat. Whenever this distinction 
is not made—i.e., when the inlet boundary conditions 
are assumed to be fixed—special care must be exercised 
in drawing conclusions from the results. For instance, 
the imaginary solutions obtained in certain cases must 
be understood to imply the physical impossibility of the 
assumed conditions rather than a limitation to the 
amount of heat that can be added. When the change 
of the upstream conditions is taken into account, the 
physical interpretation of the results becomes clear, and 
it is found that neither does the static temperature of 
the gas decrease when heat is added at an initial Mach 
Number greater than 1/+/7, nor is the amount of heat 
that can be added limited by the attainment of M = 
1.0. It is also found that the decrease of mass flow that 
may sometimes be caused by the addition of heat is in 
no way related to the attainment of sonic velocity in the 
pipe 


ASSUMPTIONS AND NOTATION 


Assumptions 


1) The flow is frictionless and one-dimensional. 
{2) The fluid obeys the equation of state for a perfect gas. 
(3) R, cp and y are constant throughout the process of heat ad- 
dition (in all numerical applications y is assumed to have the 
value 1.400). 

(4) The tube discharges directly to the surrounding atmo- 
sphere. 


Symbols 
M = Mach Number 
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6 = static temperature 
p = static pressure 
p = density 


, ol : 
T = o( + a) = stagnation temperature 


y¥-1 y/(y¥-) 
P = pl+- <a = stagnation pressure 
V = flow velocity 
A = cross-sectional area of duct 


specific heat at constant pressure 
ratio of specific heats 


7 = 

R= gas constant 

q = amount of heat added per unit mass 

AT = q/cp = increment of stagnation temperature due to the 
addition of heat 

Ss = entropy per unit mass 

m = y me = [> —— = mass flow 
R\Vo VR\VT 

k = Pi/pe 

Subscripts 

0 = outside conditions (ambient or free stream) 

t = initial conditions inside constant-area duct 

1 = modified steady-state conditions upstream of heating 
region 


= modified steady-state conditions downstream of heating 


bo 
| 


region 


UsE oF THE Macy NuMBER Functions D, G, ano N 


The analytical treatment of one-dimensional com- 
pressible flow is considerably simplified by the use of the 
Mach Number functions D, G, and N, a plot of which 
is shown in Fig. 1. With these functions, the continuity 
equation may be written in the form 


PAD/V/T = const. (1) 
In the case of isentropic flow, this simplifies to 
AD = const. (2) 


In the case of flow through a duct of constant cross 
section, Eq. (1) becomes 


PD/V/T = const. (3) 
and the momentum equation may be expressed as 
PG = const. (4) 
so that, since D/G = N, 
N/VT = const. (5) 


When TJ is constant, as across a normal shock, Eqs. 
(3) and (5) simplify to PD = const. and N = const. 


The energy equation is in all cases 
AT = q/G» (6) 
In the case of steady flow through a duct of constant 
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Fic. 1. Mach Number functions D, G, and N for air (y = 1.400). 


cross section with heat addition between two regions 1 
and 2, Eqs. (3), (4), (5), and (6) become 


P\D,/VT, = P2D2/VT2 (7) 
PG, = P2G, (8) 
Ni/VT, = No/VT2 (9) 
T, = T,; + AT (10) 


PHENOMENA RESULTING FROM THE 


ADDITION OF HEAT 


[TRANSIENT 


Consider a gas flowing in the initial condition i— 
uniform throughout-—in a duct of constant cross sec- 
tion. Now let heat be added uniformly across a sec- 
tion 7 of the duct, and let ¢ = 0 be the time at which the 
addition of heat is suddenly started (Fig. 2). As the 
gas flows through section //, it undergoes an increase 
of specific volume and therefore, because of the condi- 
tion of continuity of mass flow, an acceleration. Hence, 
the velocity of the gas will be higher in the heated re- 
gion / than in the unheated region a, and both V, and 

7, will generally be different from the initial velocity 
V,. Such a change in flow conditions can, however, 
only be established in region a through a pressure wave 
S, traveling upstream from section /7. In regions h 
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and 6, the change of flow conditions is established 
through a pressure wave S, traveling downstream from 
section /7, ahead of the interface J between heated and 
unheated gas. The upper part of Fig. 2 shows the posi- 
tion, at a time ¢ = 7, of the interface and pressure 
waves, which are represented as discontinuity surfaces 
for the sake of simplicity. The lower part of Fig. 2 is 
a time-distance representation of the propagation of 
the waves and of the displacement of the interface. 

The time-distance plane is seen to be divided into 
five regions. The state and flow parameters in the two 
regions 7 are the initial ones and therefore are known 
Also known is the amount of heat added per unit mass 
flowing through section H. The flow velocity and any 
two of the state parameters in regions a, }, and h, and 
the velocity of propagation of S;, S., and J are the un- 
knowns to be determined. This represents a total of 
twelve unknowns. The continuity, momentum, and 
energy equations applied across S;, H, J, and S» con- 
stitute the required total of twelve equations.* 

Fig. 2 only illustrates the beginning of the transient 
The pressure waves S; and S); will continue to propa- 
gate unchanged until they reach either an end of the 
duct or a change of cross section. At such points, and 
depending on the conditions prevailing there, the waves 
will generally be partly transmitted and partly re- 
flected. A similar phenomenon will occur when these 
pressure waves reach a stationary shock, as may be the 
case when they are propagated through a diffuser; in 
this case, the shock will be displaced as a result of the 
collision. 

The reflected waves will in turn be partly reflected 
and partly transmitted not only at changes of cross 
section but also at such surfaces of state discontinuity 
as have been generated earlier in the transient (e.g., 
H, I, and reflected pressure waves). This propagation 
and interaction of pressure waves and interfaces will 
continue until a new steady state is established. Then 
conditions will be uniform in each of the two regions, 


* The absence of flow through the interface reduces the equa 
tions relating regions 4 and } to the conditions that the velocity 
ind pressure be the same in these two regions and that the veloc 
ity of J be the same as that of the flow in 4 and 6. 
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Effect of heat addition on the Mach Number upstream of the heating region when the mass flow and inlet stagnation 


temperature are kept constant. 


upstream and downstream of H. These conditions, 
which will be hereafter referred to as | and 2, respec- 
tively, will generally be different not only from condi- 
tion 7 but also from the conditions a and 6 of the be 
ginning of the transient 

The final steady state depends on the configuration 
of the system, on the boundary conditions, and on the 
amount of heat added. Thus, for instance, it is found 
that, if the constant-area duct discharges directly to 
the surrounding atmosphere, the static pressure /» is 
equal to the outside pressure p, if 1/2 < 1.0 and may be 
equal to, or greater than, the outside pressure if \/. = 
1.0. Indeed, if MJ. < 1.0, any pressure difference at 
the exit would be propagated into the tube, and there- 
fore a steady state could not exist with p2 ¥ p,. If 
M, = 1.0, a steady state is possible only if pz > p,: 
the condition p, < p, would give rise to a compression 
wave, which would travel upstream into the tube since 
its velocity relative to the gas would be greater than the 
local velocity of sound;'® whereas expansion waves, 
such as would be produced by a pressure drop at the 
exit (Pp, > p,) would only travel at the local velocity of 
sound and therefore could not be propagated into the 
tube when M, = 1.0. In any case it is known that MJ, 
cannot exceed 1.0 when heat is only added to, not re- 
moved from, an initially subsonic flow.® 

If the heat addition is accomplished by combustion 
of fuel in air, the process illustrated in Fig. 2 corre- 
sponds to the idealized case of a “‘cold’’ ignition source 
distributed uniformly over a dragless flame holder 


extending over the whole cross section //, the ignition 
source being in operation before fuel is added to the air 
stream. In the case in which operation of the ignition 
source is started after the tube is already filled with 
combustible mixture, a second flame front is established 
traveling downstream into the unburned mixture; then 
the time-distance plane is divided into six regions in- 
stead of five, and the number of unknowns increases 
from 12 to 16, the four additional unknowns being the 
flow velocity and two state parameters in the added 
region and the velocity of the second flame front. The 
latter is, however, equal to the flow velocity in the re- 
gion ahead of it plus a burning velocity that is assumed 
to be known, and three additional equations (of con- 
tinuity, momentum, and energy) may be written to 
relate conditions across the additional flame front. 

A substantially similar process takes place when heat 
is added over a length of the tube rather than at one 
section, except for a gradual build-up of the pressure 
waves. 

In any case, the study of the transient shows that 
the addition of heat to a gas flowing in a tube causes 
modifications of flow and state conditions both down- 
stream and upstream of the region of heat addition. 

The main usefulness of a quantitative analysis of the 
transient phenomena described above would lie in the 
information that it would provideontherate of approach 
to the final modified steady state, as affected by size, 
configuration, and boundary conditions. It is possible, 
however, to arrive at a determination of the modified 
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steady state directly—i.e., without consideration of the 
transient—and this will be done in the following sec- 
tions. 


THE MOpIFIED STEADY STATE 


To illustrate the manner in which the modified 
steady state can be determined directly from the known 
initial and boundary conditions and from the amount of 
heat added, two simple cases will be treated here: 
namely, the cases in which either the mass flow and 
stagnation temperature at the inlet or the stagnation 
temperature and stagnation pressure at the inlet re- 
main constant as heat isadded. Both cases are of prac- 
tical importance: the first one occurs, e.g., when the 
air enters the pipe through a diffuser with a supersonic 
inlet; the second one occurs when the air is fed into 
the pipe from a constant pressure reservoir through a 
diffuser in which the flow is everywhere subsonic. 

The treatment will be limited here to the determina- 
tion of Mf, and J/; as functions of the initial conditions 
and of the amount of heat added. All other parameters 
pertaining to the modified steady state may then be 


1949 


computed directly from the equations defining these 
parameters and from Eqs. (7) and (8). 


(a) Constant m and T, 


Since 7; = 7;, the condition that the mass flow be 
the same in condition 1 as in condition 7 may be written 
in the form 


P,\D, = P,D, = const. (11) 
From Eqs. (7) and (11) one obtains 
P.Di/VT, = P2D2/ VT, 


i.e., 


pM, y: + 1Stue/ VT. = 
pile + 1 as | VT; 


If Mz < 1.0, then pz = pp, and, since this investigation 
is limited to subsonic initial conditions, £; = p,, hence, 


pb, = py Then, solving for A/s, one obtains 


ais a a 
M,? = — fyi + + 27 — pate 4 + tly Ni + 4) - | (12 
yo T; 


Eq. (12) shows that M], increases continuously with 
AT. When Eq. (12) gives Mz > 1.0, the condition 
P2 = p, does not hold any longer, but then Eq. (12) is 
replaced by the condition 1/, = 1.0. 

The determination of 1/4, may then be made from 
Eqs. (9) and (12). 

The values of the Mach Numbers upstream and 
downstream of the region of heat addition in the modi- 
fied steady state, as computed from Eqs. (9) and (12), 
are plotted in Figs. 3 and 4 for various initial Mach 
Numbers and amounts of heat added. It is seen that 
the upstream conditions are considerably modified by 
the addition of heat. When sonic velocity is reached 
at the exit (1/2 = 1.0), the Mach Number at the inlet 
continues to decrease as more heat is added but in a 
manner that is then independent of the initial Mach 
Number (top curve of Fig. 3). The corresponding 
changes of P; may be calculated from Eq. (11), which 
shows that, since 4, decreases, P; will increase as heat 
is added. 


(b) Constant P, and T; 


P, = P, = kp, and Eq. (8) can be written in the 
form 


kpoGi = pol + yM,?) 
whence, if M/; < 1.0 and therefore p. = p,, 


Mz = V(1/y)(kG — 1) (13) 


Eq. (13) 


When the exit velocity is sonic, is replaced by 


the condition MJ, = 1.0. 

















Effect of heat addition on the Mach Number down- 
flow and inlet 


Fic. 4. 
stream of the heating region when the mass 
stagnation temperature are kept constant. 
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Fic. 5. Effect of heat addition on the Mach Number up- 
stream of the heating region when the stagnation pressure and 
temperature at the inlet are kept constant. 


From Eqs. (9) and (10) one obtains 


AT/T, = AT/T, = (N2/N,)? - (14) 
) and (14) relate 14, and ., to the amount 
of heat added. These relations are shown in graphical 
form in Figs. 5 and 6. The effect of heat addition on 
upstream conditions is even more marked in this case 
than in case (a). Of particular interest is the fact that, 
unless k > [(y + 1)/2]”~-” 1.893, the down- 
stream Mach Number (\Me) always decreases as heat is 


Eqs. (13 





° + 4 6 8 
AT/Tj 


Fic. 6. Effect of heat addition on the Mach Number down- 
Stream of the heating region when the stagnation pressure and 
temperature at the inlet are kept constant. 


GAS FLOWING IN A PIPE 89 


added and id approaches the asymptotic value M,; = 
JV (k — 1)/y; this can be verified from inspection of 
Eqs. (13) and (14) and noting that, as AT/T,; > @, 
N, — 0 and therefore G; — 1.0. 

If k > 1.893, M2 remains constant and equal to 1.0 
until a certain amount of heat has been added, this 
amount depending on the value of &. Further addi- 
tion of heat will again result in a decrease of M2. Only 
ifk > y + 1 = 2.4 does M2 remain constant and equal 
to 1.0, regardless of the amount of heat that is added. 
The value k = 1.893 corresponds to the critical pres- 
sure ratio that would make M/, = 1.0. 

The fact that MM, decreases when heat is added must 
not be interpreted to mean that the gas velocity V2 
also decreases. From Eq. (14) and from the relation 
V2 « M2v/6, one obtains 

MV + (AT/T)_ M? 


f —=—=—<—<———===== ry ————- (15) 
“V1 + Uy — 1/2]M My ~ Nil + 7M") 
Logarithmic differentiation gives 
dV2 2dM2 dN, (16) 
Ve ~ M(1 + M2? :) Ni 








The second ratio on the right-hand side may be written 
as 


dN, 


dN, on ee) 
7 (+) () (= am,) 2: (16a) 


If Mz < 1.0, then, using the definition of the Mach 
Number functions G and N and Eqs. (13) and (16), 
dV2 2 
——~ 2 2 (17) 
Eq. (14) shows that N2 > N, and therefore M, > M,. 
On the other hand, it has been shown above that, when 
If M, = 1.0, Eq. (15) shows directly that V2 in- 
Therefore one may conclude that the gas velocity in 
region 2 always increases when heat is added if P, and 
in Fig. 7. Under the same conditions it is seen from 
the definition of mass flow that the latter must de- 


and (16a), one obtains 
E (+) |= 
V2 kG, My 2 Ve 
M. < 1.0, dM2/Msz < 0. Therefore dV2/V2 > 0. 
creases as Al + (AT/T)). 
T, are constant. This effect is graphically illustrated 
crease, since //, decreases as heat is added. 


DISCUSSION 


The determination of the new steady state in any 
specific case can always be made by a procedure similar 
to that illustrated in the preceding section—i.e., from 
the amount of heat added, the initial and boundary con- 
ditions (which include the characteristics of the air 
supply system, described by two equations relating 
m, T, and P;), and from Eqs. (7) to (10). 

The boundary conditions—in particular the charac- 
teristics of the air supply—may vary so widely from 
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Fic. 7. Effect of heat addition on the exit velocity when the 
stagnation pressure and temperature at the inlet are kept constant. 


system to system as to make it impossible to determine 
the modified steady state in a general form. There 
are, however, certain aspects of this problem that can be 
treated—at least in a qualitative manner—independ- 
ently of the boundary conditions. A general treatment 
of this sort will be applied here to a discussion of some 
of the conclusions that have been reached by a number 
of previous investigators—-namely, that: 

(a) When heat is added to a gas flowing at an initial 
Mach Number within the range 1/Wy < M < 1.0, 
the static temperature of the gas decreases. 

(b) When heat is added to a subsonic flow, the Mach 
Number always increases, approaching the sonic state. 

(c) Heat cannot be added to a gas if its initial veloc- 
ity is equal to the velocity of sound, and therefore the 
condition 1J = 1 sets a definite limit on the amount of 
heat that can be added. 

(d) When 1/7 = 1.0 is attained, any further addition 
of heat will result in a decrease of mass flow. 

These conclusions will be reviewed in the following 


sections. 
Effect of Heat Addition on Static Temperature 


The statement has been made that a decrease of 
static temperature, when heat is added at an initial 
Mach Number within the range 1/-~/y < M < 1.0, 
would not be in contradiction to the second law of 
Thermodynamics. That such a decrease would ac- 
tually be a violation of the second law can, however, 
be shown as follows. 

The entropy of the gas in region 2 may be measured 
using as level of reference that of an upstream region 
where no change of entropy will result from the addi- 
tion of heat in the tube. This condition will generally 
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be found to be satisfied a short distance upstream of the 
entrance to the air supply system. Since the entropy 
level of the medium in the chosen region of reference 
is the same regardless of the amount of heat that is 
added in the tube, it follows that the entropy per unit 
mass in region 2 is determined by any two of the state 
parameters in this region. Thus, 


Sg —- 55 = & log (02, es R log (p2/ Po) 


where so and (c, log 6, — R log p,) are independent of 
the amount of heat added in the tube. 

In all cases of conceivable interest, d72/d( AT) > 0— 
i.e., the amount of heat added in the tube is not ex- 
ceeded by a simultaneous heat removal in the intake 
system. Then, an increase of AT can only be achieved 
by an increase of the net amount of energy added per 
unit mass flow of the gas between regions 0 and 2, and 
therefore the entropy level in region 2 must increase as 
AT increases. 

If \J, < 1.0 and therefore p2 = const., 

dsz/d( AT) = (¢,/02)(d@2./d(AT)] > 0 


whence d@,/d( AT) > 0. 

If Mz = 1.0, then 7, = [(y + 1)/2], and the con- 
dition d72/d(AT) > 0 again yields d@./d(AT) > 0. 

Therefore, except in the irrelevant case in which the 
addition of heat is more than completely offset by a 
simultaneous heat removal, the addition of heat will al- 
ways result in an increase of static temperature of the gas 
above its initial value, regardless of the initial Mach 
Number. 

Three illustrations of this effect are given in Figs. 8, 
9, and 10, which refer to specific examples. The mass 
flow and stagnation temperature at the inlet are con- 
stant, and the initial Mach Numbers are 0.4, 0.9, and 
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Fic. 8. Effect of heat addition on static temperature when 
the initial Mach Number is 0.4 (mass flow and inlet stagnation 
temperature constant). 
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Fic. 9. Effect of heat addition on static temperature when the 
initial Mach Number is 0.9 (mass flow and inlet stagnation 
temperature constant). 


1.0, respectively. Curves of static temperature versus 
Mach Number within the heating region and for various 
steady régimes of heat addition are drawn in all three 
The curves 1 and 2 relate static temperatures 
and Mach Numbers in regions 1 and 2. Before heat is 
added, the condition is that indicated by point 7. Re- 
ferring to Fig. 8, when an amount of heat corresponding 
to, say, AT/T, = 3 is added to the gas, condition 1 be- 


figures. 


comes that corresponding to point A and condition 2 
between—i.e., 
Number and 


that corresponding to point B. In 
within the heating region—the Mach 
static temperature are related as shown by the curve 
labeled AT/T,; = 3. It will be noted that the Mach 
Number in region 1 has decreased from 0.4 to 0.24, 
whereas the Mach Number in region 2 has increased 
from 0.4 to 0.77. The static temperature has increased 
only slightly in region 1 but considerably in region 2, 
As more heat is added, the static temperature 62 con- 
tinues to increase, not only up to Mz = 1/7/¥ but also 
beyond this point and after 1/2 = 1.0 has been reached. 
For instance, when a further increment of heat from 
AT/T, = 4 to AT/T; = 6 is added after MJ has 
teached the value 1/+/y (point C), the condition in 
tegion 2 changes from C to D, indicating an increase of 
static temperature. An analysis based on the assump- 
tion of fixed inlet conditions would lead to the belief 
that the change of condition 2 would take place in this 
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case from E to D and that therefore @: would decrease. 
This appears to be the basic reason for the already men- 
tioned misinterpretation of the effect of heat addition 
on static temperature. It is seen from Fig. 8 that no 
significant change in the effect of heat addition occurs 
either at Mz = 1/~/y or at Mz = 1.0. 

The relation between @ and .\/ within the region of heat 
addition for various values of A7’/7,; has been derived 
by various authors and may be written in the form 
[obtainable from Eq. (5)] 


Ni/-VT, = M ‘\/o(1 + yM?) = const. 


For any given J/,, values of J/, corresponding to 
various values of A7T/7T, may be taken from Fig. 3. 
Differentiation of Eq. (18) with respect to 1/ shows 
that @ is maximum where MW = 1//y. Hence, 6; may 
be greater than #. This may be interpreted in the 
following manner: The gas, as it flows in steady ré- 
gime through the region of heat addition, undergoes a 
pressure drop—i.e., an éxpansion. Therefore, the 
static temperature rise of the air as it crosses the heat- 
ing region is lower than it would be if heating were 
taking place at constant pressure and may even turn 
In other words, the increase 


(18) 


into a temperature drop. 
of static temperature resulting from the addition of 
heat may be greater in region 1 than in region 2. Eq. 
(9), written in the form 


M,/V0(1 + yh?) = 


M2/V00(1 + yM2*) 


6/T; 


o9 
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Fic. 10. Effect of heat addition on static temperature when 


the initial Mach Number is 1.0 (mass flow and inlet stagnation 
temperature constant). 
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shows that 06,/62 > 1.0 if 
(1 + yM,*)/Mz > (1 + yM,2)/M, 


It can easily be verified that, since M; > M, [Egs. 
(9) and (10)], this condition may occur if M, > 1/y 
and will always occur when Mf, > 1/+7/y. 

Two such cases are shown in Figs. 9 and 10 where it 
is seen that for extremely small amounts of heat added 
the static temperature in region 1 may become higher 
than in region 2. 

It is well known® that, whereas the static tempera- 
ture may increase and decrease across the region of 
heat addition in a steady régime as the Mach Number 
increases, the entropy increases continuously and only 
reaches a maximum at M; = 1.0. This can also be 
verified directly as follows. The entropy level, s, at 
any point in the region of heat addition may be ex- 
pressed in the form 


S= 53, a Cp log te T;) — R log (P P,) 
Then, from Eqs. (+) and (5) one obtains 
s = const. + 2c, log N + RlogG 


showing that s is maximum at JJ = 1.0, since the maxi- 
mum of both N and G occurs at this Mach Number 
(Fig. 1). Therefore, the second law of Thermody- 
namics calls for @ > 0, but not for @. > 4;. 


Effect of Heat Addition on the Mach Number of the 
Heated Flow 


If M, < 1.0 and therefore p. = const., then the den- 
sity of the gas in region 2 is inversely proportional to the 
static temperature #. Then, 


m= pV, ( V2/2) 


and logarithmic differentiation with respect to AT yields 


1 dm — = 1 dV. 1 dO, 19 
md(AT) V2d(AT) 62d(AT) (19) 
From the definition of Mach Number, one obtains 
, aMy = 1 dV, a P d02 20 
M, d( AT) V2 d( AT) 26. d( AT) (20) 
Eqs. (19) and (20) give 
1 dM, 1 dm l dO. 
= ~ (2) 
M2 d( AT) m d( AT) 20, d( AT) 


Since d@,/d(AT) > 0, it follows that, if m = const., 
the Mach Number in the heated region must always 
increase as heat is added. 

If, however, the characteristics of the air supply or 
intake system are such that the mass flow decreases as 
heat is added, then M; may increase or decrease de- 
pending on the relative rate of change of m and 62. 

Two special cases of constant and of decreasing mass 
flow have been treated in detail in the discussion of the 
modified steady state, and it has already been pointed 
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out that the addition of heat would cause an increase 
of M, in the former case and a decrease in the latter. 

Therefore, the Mach Number of a heated subsonic 
flow may increase or decrease as heat ts added, depending 
on the boundary conditions. 


Addition of Heat to a Sonic Flow 


The case in which the initial Mach Number (1/,) is 
1.0 will be discussed first. This is the case in which, 
according to some of the previous investigators, no ad- 
dition of heat would be possible. 

Eq. (9) shows that 1/; must always be lower than 
M.. Since M2 cannot exceed 1.0, it follows that, if 
M, = 1.0, the transient resulting from the addition of 
heat will modify conditions upstream of the region of 
heat addition so that, in the modified steady state, 
M, < 1.0. In other words, heat can be added to a flow 
that is initially sonic, but M, is never equal to 1.0. 

Failure to distinguish between .1/; and .l/, has ap. 
parently been the cause for the widespread belief that 
heat cannot be added to an originally sonic flow. 

A similar argument may be applied to the case in 
which heat is added to a flow that was initially sub- 
sonic but has already been heated until 1/2 has reached 
1.0. In this case, it follows from the conclusion of the 
preceding paragraph that MM, must already be lower 
than 1.0 even before the further addition of heat. The 
effect of this additional heating will merely be to con- 
tinue to change the conditions in region 1. The Mach 
Number in region 2 may remain constant or even de- 
crease as heat is added, again depending on the bound- 
ary conditions. 

The effect of heat addition to an initially subsonic 
flow is illustrated in Figs. 8 and 9. The case of an 
initially sonic flow is shown in Fig. 10. All three ex- 
amples refer to configurations with intake systems that 
maintain a constant mass flow. 

The assumption of a constant mass flow is not in- 
consistent with the condition 1/7, = 1.0 or JJ, = 1.0. 
The mass flow depends entirely on the conditions up- 
stream of the region of heat addition and on the charac- 
teristics of the intake system. The effect of heat addi- 
tion on conditions in region 1 is substantially the same 
whether the heat is added when J/, = 1.0 or when J/2 < 
1.0 (see, e.g., Figs. 3 and 5 or 8,9,and 10). The often- 
made statement that the mass flow must decrease 
when the flow becomes sonic in region 2 appears to be 
the result of an attempt to interpret one of the singu- 
larities to which the theory inevitably leads at J/ = 
1.0 whenever the inlet conditions are assumed to be 
fixed. On the other hand, if the air intake system is 
one in which any increase of back pressure is always 
accompanied by a decrease of mass flow, then this de- 
crease will occur as soon as heat is added, regardless of 
the value of the Mach Number in region 2. This 
shows that the changes in mass flow that may accompany 
the addition of heat are in no way related to the attainment 
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Example 1. 


of sonic velocity in the tube. The examples in the next 
section will illustrate this fact. 


EXAMPLES 


The following examples will illustrate the manner in 
which the procedures described in the preceding sec- 
tions can be used to calculate the state and flow modi- 
fications caused by heat addition. In all three ex- 
amples the air enters the pipe through a ram intake, 
consisting of a simple diffuser. Therefore, the stagna- 
tion temperature 7) is constant in each case and equal 
to the initial stagnation temperature 7, which in turn 
is equal to the free-stream stagnation temperature 7,. 


Example 1 


The free-stream Mach Number is M, = 1.5, and the 
air enters through a stable supersonic diffuser®® de- 
signed for this Mach Number. The configuration is 
shown in Fig. 11. In the absence of heat addition, the 
initial Mach Number in the pipe (1/,) and the location 
of the shock in the diffuser are determined as follows. 
Since the flow is subsonic in the pipe, p; = ~,. Then 
the continuity equation [Eq. (1)], applied between the 
intake and the pipe (with 7, = 7,) and solved for 
M,, yields M, = 0.4. The ratio P,/P, may then be 
calculated, since p; = p, and both M, and MM, are 
known. Since the flow is assumed to be frictionless, 
the loss of stagnation pressure may be entirely at- 
tributed to the shock, and the Mach Number at which 
this shock occurs may be obtained from the Rankine- 
Hugoniot equations or from tables.*! Knowing the 
Mach Number at which the shock occurs, the cross- 
sectional area of the diffuser at the section at which the 
shock is located may be calculated from Eq. (2), ap- 
plied between this section and the intake. The loca- 
tion of the shock (x) is thus determined for the initial 
condition (AT/T, = 0). When heat is added, the 
values of the Mach Numbers M,; and M2 may be ob- 


GAS FLOWING IN A PIPE 93 


tained from Figs. 3 and 4 or calculated as explained in 
the discussion of the modified steady state for the case 
of constant mass flow. The manner in which J, 
changes as heat is added is shown in Fig. 11. As 
long as .f2 < 1.0, the shock location (x) may be deter- 
mined from Eq. (8), which gives P; and therefore the 
pressure ratio P;/P,. As the amount of heat added in- 
creases, 4, increases until it reaches 1.0. At this 
point (A7T/T, = 6.3), the shock is still well inside the 
diverging part of the diffuser, as can be seen from curve 
x of the figure, the ordinates of which can be measured 
directly on the sketch at right. Further addition of 
heat results in further displacement of the shock without 
change in A/,. The shock location for various values 
of AT/T, > 6.3 (i.e., with Mz. = 1.0) can be deter- 
mined from the known values of J/, (Fig. 3) and from 
Eq. (2), applied between the intake and the section at 
which the shock is located and then between the latter 
and region 1, together with the condition that PD is 
constant across the shock: It is seen that, up to the 
high values of stagnation temperature rise that are 
considered in this example—as high as or higher than 
may be obtained from the combustion of ordinary fuels 
—the shock remains within the diffuser and therefore 
the mass flow remains constant. This example shows 
that the attainment of sonic velocity at the exit of the 
pipe does not introduce any discontinuity in the effect 
of heat addition, except for the rate of change of M2. 


Example 2 


The free-stream Mach Number is 0.8, and the air is 
taken in through a simply diverging diffuser (Fig. 12). 
In the absence of friction losses, P; = P, = P,, and the 
initial Mach Number in the pipe is equal to the free 
Mach Number (M, = 0.8). Since P; and 7; are con- 
stant, the values of 4/, and JJ, for various amounts of 
heat added can be obtained directly from Figs. 5 and 
6 or calculated as explained in the discussion of the 
modified steady state for this case. Again, since P; 
and 7; are constant, the mass flow is simply propor- 
tional to D, (see definition of m) and can therefore be 
computed from the known values of 1/;. The manner 
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in which the mass flow and the exit Mach Number vary 
as heat is added is shown in Fig. 12. This is clearly a 
case in which the Mach Number at the exit not only 
does not approach unity but actually decreases when 
heat isadded. Also, the mass flow is seen to be decreas- 
ing continuously as heat is added, even though the 
flow is never sonic anywhere in the pipe. 


Example 3 


The free-stream Mach Number is the same as for 
Example 2 (.\/, = 0.8), but the diffuser is so designed 
that the flow becomes supersonic over a certain region 
downstream of the throat (Fig. 13). The manner in 
which the shock location, 1/2, and the mass flow change 
as heat is added is shown in Fig. 13. The initial Mach 
Number .)/, and shock position are determined as in 
Example 1, as are \/; and .W/, and the shock location 
corresponding to various amounts of heat added. The 
mass flow remains constant until the shock reaches the 
throat—i.e., as long as A7’/7; < 1.84. In this range, 
My increases from 0.4 to 0.654. Beyond this range, 
the flow is everywhere subsonic, and the computation 
proceeds as in Example 2. Then both the mass flow 
and the Mach Number decrease. Again, in this case 
it is seen that the appearance of discontinuities in the 
rate of change of flow conditions is in no way related to 
the attainment of sonic velocity in the pipe. 


CONCLUSIONS 


The following conclusions may be drawn on the 
changes of state and flow parameters that are caused 
by the addition of heat to a subsonic or sonic flow as 
considered in the preceding analysis: 

(1) The addition of heat has the effect of modifying 
the state and flow conditions not only downstream but 
also upstream of the region of heat addition. These 
modifications are often greater upstream than down- 
stream of the heating region. The nature of these 
changes depends on the boundary conditions. 

(2) The Mach Number of the flow downstream of the 
heating region does not necessarily increase when heat is 
added. Instead, it may increase or decrease, depend- 
ing on the boundary conditions. 

(3) The effect that the addition of heat may have on 
the mass flow through the pipe depends exclusively on 
the boundary conditions and is in no way related to the 
attainmeént of sonic velocity in the heated region. 

(4) There is no aerothermodynamic limitation to 
the amount of heat that can be added to a gas flowing 
in a pipe. The attainment of sonic velocity in the 
heated region does not limit this amount. Also, it is 
always possible to add heat to a gas that is initially 
flowing at sonic velocity. However, the Mach Num- 
ber upstream of the heating region will always become 
lower than 1.0 as a result of the addition of heat. 

(5) The static temperature of a flowing gas does not 
decrease when heat is added to it at an initial Mach 
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Number within the range 1/+ “'y < VW, < 1.0; instead, 
it always increases, although the temperature rise due 
to heat addition may, in this range of Mach Numbers, 
be greater upstream than downstream of the heating 


region, 
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Applications of the Theory of Free Molecule 
Flow to Aeronautics 


HOLT ASHLEY* 
Massachusetts Institute of Technology 


SUMMARY 


The breakdown of certain assumptions of conventional aero- 
dynamics in the case of highly rarefied gas is shown to call for the 
use of the theory of free molecule flow in studying forces on bodies 
moving above 120-150 km. in the earth’s atmosphere. A tenta- 
tive standard table of atmospheric properties from 120-220 km. 
iscomputed. The parameter M o, ratio of flight speed to most 
probable speed of air molecules, becomes useful in comparing free 
molecule flows. 

Techniques are developed for finding analytically or numer- 
ically the lift and drag coefficients of any object, based on two 
extreme opposite assumptions concerning the nature of reflection 


of molecules from surfaces. These methods are employed to find 


coefficients of a number of simple bodies. 

The drag of missiles is demonstrated to be negligibly small 
ibove 120 km. except at extremely high speeds. The possibility 
of aerodynamically sustained flight is discussed. 


INTRODUCTION 


A’ A RESULT OF DEVELOPMENT of propulsive devices 
involving the flow of fluid at very high densities 
and temperatures and of the penetration of missiles to 
extreme altitudes in the atmosphere, certain new as- 
sumptions about the nature of gases must sometimes 
be introduced into aerodynamic theory. If compress- 
ible gas flows are considered, the ideas of a continuous 
medium and a perfect gas with constant y = C,/C, 
will be found to underly most developments. These 
assumptions are valid for almost every problem treated 
in ‘classical’ aerodynamics; errors introduced by them 
are small compared with inaccuracies due to the postu- 
lation of inviscidity, linearization of differential equa- 
tions, etc. 

Let us reconsider the mathematical “‘model”’ for very 
high and low-density cases. In high-density and tem- 
perature flow there is no objection to the continuous 
medium idea, but the state of the gas must be repre- 
sented by one of those equations developed to describe 
affairs as the critical point is approached. Very low- 
density aeronautical problems were first treated by 
Zahm! in 1934. He gave the title “‘superaerodynamics’”’ 
to the field of mechanics of highly rarefied gases. 
Tsien? has delineated the field more precisely by estab- 
lishing regions in terms of two familiar parameters of 

Complete paper received February 12, 1948. Condensed ver- 
sion received October 11, 1948. Complete report was published 
as a Sherman M. Fairchild Publication Fund paper and can be 
ordered as Preprint No. 164 from the Preprint Department, I.A.S., 
2 East 64th St., New York 21, N.Y. 

* Assistant Professor of Aeronautical Engineering. 


aeronautical engineering and describing what model is 
to be used to represent the gas in each of these regions. 

For extremely rarefied gas, corresponding roughly to 
altitudes above 150 km., superaerodynamics employs 
the idea of free molecule flow; no longer can the gas 
be justifiably termed continuous, and flows are de- 
scribed by powerful and experimentally well-verified 
principles of the kinetic theory of gases. The present 
paper is devoted to the study of flows over bodies in 
this so-called ‘‘free molecule region.’’ Techniques for 
computation of forces and moments on simple geo- 
metrical objects and on bodies of arbitrary shape are 
outlined. Actual forces given were computed using ap- 
proximate mean properties of the high atmosphere 
listed in Table 1. These were computed by the 
N.A.C.A. method’ from the assumed temperatures and 
compositions indicated. 

Tsien defines the region of applicability of free mole 
cule theory as corresponding to 


\/L > 10 
where 


AX = mean free path of gas molecules 
L = typical length of object in flow 


Thus only altitude and body dimensions count in as- 
signing a case to this range, and for ordinary objects 
of 10° cm. size all flight above 170,000 m. in the atmos- 
phere should be strictly so assigned. In the hope that 
experimental results may later justify extrapolation to 
somewhat lower altitudes, this paper considers heights 
between 120,000 and 220,000 m. The extreme thinness 
of the atmosphere above 220,000 m. suggests that it 
will have negligible influence on any flight there. 
Flight at 220,000 m. at the fantastic speed of 10,000 m. 
per sec. corresponds to a dynamic pressure of only 
about 40 dynes per sq.cm., while static pressure at the 
earth’s surface is 10° dynes per sq.cm. 

Two final problems must be considered before treat- 
ing free molecule flows in detail. The first involves 
the nature of the distribution of velocity of thermal agi- 
tation among the molecules of the gas. The so-called 
Maxwell or normal distribution can often be presumed 
with great precision, but it is maintained by just that 
mechanism of molecular collisions which is considered 
negligible in the free molecule region. However, the 
absence of agencies to destroy the distribution, its 
maintenance by the detailed balancing operation, and 
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— a © Te 2 TABLE 1 fs oe £ 
Estimated Mean Properties of the Atmosphere Above 120 Km. Altitude H 
ee ny, € 
Atmosphere Consisting of Atomic Oxygen (O) and Molecular Nitrogen (N2) is | 
Altitude Pressure, p Density p X N 
(M.) (Dynes per Sq.Cm.) (Gm. per Cu.Cm.) (Cm.) {Molecules per Cu.Cm.,) 
120,000 0.4707 36.15 X 10-4 19.90 90.88 X 10% dN 
130,000 0.2322 16.99 42.21 42.72 
140,000 0.1185 8.279 86.35 20.82 
150,000 0.06241 4.175 170.7 10.50 
160,000 0.03383 2.170 327.5 5.455 Te 
170,000 0.01884 1.161 610.3 2.918 c 
180,000 0.01075 0.6369 1,109 1.601 the 
190,000 0.006268 0.3578 1,968 0.8996 ; 
200,000 0.003737 0.2058 3,411 0.5173 - 
210,000 0.002272 0.1208 5,792 0.3038 
220,000 0.001405 0.07224 9,657 0.1816 
Atmosphere Consisting of Atomic O and N 
Temperature are 
Altitude Pressure p Density p r N (Degrees) ha 
(M.) (Dynes per Sq.Cm.) (Gm. per Cu.Cm.) (Cm.) (Molecules per Cu.Cm.) Kelvin ‘ 
120,000 0.4707 36.15 XK 107" 19.90 150.3 X 10-# 375.0 TI 
130,000 0.3070 13.59 31.92 93.40 393.5 m¢ 
140,000 0.2044 8.638 50.05 59.38 412.0 
150,000 0. 1387 5.611 76.82 38.58 430.5 a 
160,000 0.09578 3.715 115.7 25. 54 449.0 gre 
170,000 0.06723, 2.504 yg ie 17.22 467.5 ; 
180,000 0.04788 1.716 248.9 11.79 486.0 = 
190,000 0.03456 1.193 356.9 8.200 504.5 
200,000 0.02527 0.8415 504.4 5.785 523.0 
210,000 0.01871 0.6016 703.5 4.136 541.5 " 
220,000 0.01399 0.4349 970.0 2.990 560.0 
the continual upward migration of normally distributed tabular integration of force distribution will always pro- . 
molecules from below, all urge the unqualified accept- duce the desired coefficient. 
ance of Maxwell's result. The method is first demonstrated by reviewing the 
A second and much more difficult problem concerns — shear and normal forces acting on unit area of a flat sur- 
the mechanism of reflection of molecules incident upon _ face oriented arbitrarily in a stream of velocity U. This , 
a body in the flow. This is the critical question in fix- problem is solved most concisely by Tsien? and is also 
ing the magnitude of forces and the efficiency of lifting treated with the introduction of several new functions 
surfaces, and it has not yet been answered experiment- by Snow.‘ Normal and shear stresses are separated 
ally with finality. It will be discussed fully below. into those due to the incoming stream of molecules (p, - 
Associated with the reflection problem is that of energy and r;) and those due to the reflected stream (p, and 
balance between surface and surroundings. t,). The former are developed, as follows: Te 
Let &’, n’, and ¢’ be velocity components of a molecule 
METHOD OF COMPUTING FoRCE AND MoMENT in directions x’, y,’ and 2’ fixed relative to the mean I 
° (¢ 
COEFFICIENTS OF AN OBJECT IN FREE MOLECULE motion of the gas. Impacts between molecules and 
FLow mutual forces are neglected, so the dynamic theory of 
7 gas states that the distribution of velocities will be 
_ The computation of actual forces and moments act- Maxwellian—that is, the number of molecules per unit 
ing on a body moving in the free molecule region of the volume having velocities in the range ¢’ > £’ + d¥’, 
atmosphere requires knowledge of the state of the air 4’ —, n’ + dn’, ¢’ > ¢ + dt’ is 
at flight altitude and the shape of the body. For A\*/ no 
we ; coe i _- , mn 
many simple geometrical shapes, force coefficients can dN’ =N 5 eT ME + +E) Gerdnide’ — (1) pr 
be found analytically in terms of ‘‘free molecule Mach me 0. 
Number” ./,, and certain other parameters. Atten- where 
tion in this paper is mainly focused on finding lift and N = number of molecules per unit volume 4 
drag coefficients C; and Cp, defined as usual and referred h = 1/2RT, = 1/c2 
to a typical area A, and dynamic pressure (!/2)pU?. ; ’ ; 
Consider unit area of a plane surface moving with 
M. = U/C : T T: . : a ya 
velocity components e, U, e2.U, e;U in directions x’, y’, 2’. 
where C = V2RT is the most probable velocity of (a? +e? +e? = 1) W 
the molecules. mn 7 
By obvious expansion of the technique of this sec- Fix coordinate system x, y, 2 in the surface parallel to 
tion, any other force or moment coefficient may be x’, y’, and 2’ and such that x is the outward normal to wl 
computed. For objects of complicated shape or of the surface. Then any molecule has velocity compo- Ch 
shape inexpressible with elementary functions, a simple nents in x, y, 2 given by fre 
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f= —-aU, n=9'-eU, $= $'—eU 


Hence, the distribution of molecular velocities in x, y, 2 
is given by 


Kd 
dN =WN (*) ‘eae + e:U)? + (n + e2U)? + (F + e3U)?] x 
T 


dtdndt (2) 


' To find the number of molecules n striking unit area of 


the plane surface in unit time, note that such molecules 
in velocity range 

f7o+ df 

are contained in a cylinder based on the surface and 
having slant height VE + n* + ¢ and altitude —é. 
The minus sign multiplies the altitude because only 
molecules with x velocities in the negative range 
—«o < §< Ocan strike the plate. wm is then the inte- 
grated sum of —£ dN over all possible velocities of 


gé—> &+ dé, n> + dn, 





impinging molecules. 


*) 
n= -v(“) of af anf” x 


ele + e1U)? + (n + e2U)? + (F + e3U)?] tdt 


—eUth eU r , \ 
— — [i 
M5 oo - [1 + erf (e:U+/h)] 
(3) 


“erf’’ represents the conventional error function of sta- 
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tistical theory, defined by 
erf (w) = e/va) fr e~* ds 
0 


This is an odd function and increases from 0 when w = 
0 to 1.0000 when w > 3; a good table may be found 
in Jahnke and Emde.® 

If the plane surface moves at an angle of attack a to 
the stream, 


é€} = sin a, @& = COS a, é3 = 0 
Furthermore, 
UVh = U/V2RT, = U/C, = Me 


where subscript 7 refers to the impinging stream. Eq. 
(3) may then be written 


NC, { a 
2 Vr 





n = — Mo? sin’ a Ae 


M,. sin a[l + erf (MV, sin ol} 


If molecules have average mass m, then mN = p, and 
the mass m, of stream striking unit area in unit time is 


ms 1 = 
mM, = “= es ¢ 


Tv 


M ~? sin? @ 4 


M,. sin a [1 + erf (7, sin a)}} (4) 


To find the forces on the surface, consider first the total momentum of the stream m, in some arbitrary direction 
specified by e:’, @2’, es’. This is found by integrating the momenta of molecules in each velocity range as follows: 


(or’, oe") ~*~ 


Ties’, ex’, e:’) = et 


2hu 


The entire momentum of the stream in a direction 
normal to the surface is transferred to the surface as 
pressure p;. It is found by setting e;’ = —1, @’ = e;’ = 
0. If, as before, e; = sin a and e, = cos a, Eq. (5) gives 


bs sin a 


(1/,)pU? Max 


1 , 
(sr + sin? «fh + erf (M,, sin a)] (6) 


—Mo!? sin’ a of. 





0, this result reduces to 
(1/2)pRT, 


which is just half the total pressure of a gas at rest anti- 
cipated by tue perfect gas law. The other half derives 
from elastic reflection of the molecules. 


When U = 


= ('/,)pc? = 


E ro + €1(€1€1’ + €2€2’ + ew) | + erf euvi) | | 


eed f- anf E(e,’ + €2'n + es’ ¢)e~"E + e:U)? + (9 + e2U)? + (F + e2U)?] dé 


7, (eres + €2€2' + e3€3 lata + (5) 


| 


/ 





If the stream of gas is first adsorbed by the surface 
before reflection, its entire tangential momentum is 
also transferred. The shear stress 7, thus produced 
may be derived from Eq. (5) by setting e’ = —1, 4’ = | 
3" = (0: 

T; COs @ 


a —Mo! sin? a 
(pl? = MeV + 


sin a cos a [1 + erf (M, sin a)}] (7) 


Additional stresses on the plane surface due to reflec- 
tion of the molecules depend on the mechanism of re- 
flection. This paper will consider two types of reflec- 
tion; they may be considered as extremes between 
which the actual process lies. The concepts of specular 
and diffuse reflection have been recognized for many 
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years, since the studies of highly rarefied gases by hard to imagine much truly elastic reflection of small The 
Knudsen and others. They are discussed at length in gas molecules from that surface. to 0 
relation to the present problem by Sanger.® In the (b) Diffuse Reflection—The molecules are trapped 
two cases, stresses are found as follows: into the interstices of the surface and re-emitted dif- 
(a) Specular or Elastic Reflection—vThe molecules fusely in directions independent of those of incidence, 
are assumed to rebound perfectly elastically from the The temperature of the molecules is adjusted toward 
surface, like light rays from a mirror. There is no that of the surface in a manner to be discussed later. 
adjustment of the gas temperature toward that of the The temperature of the re-emitted stream is designated Wh. 
surface. Since the angle and speed of reflection are 7,. There is no preferential direction of re-emission, 
equal to those of incidence, no tangential momentum sor, = 0, and the total tangential stress 
is transferred to the surface and the total shear stress whi 
ei : e ; = g 
7=0. Itisimplied that7, = —7; Since the normal ™ \ kine 
momentum of the gas stream is just reversed, the total Texts on dynamic theory of gases’ show that the 
pressure becomes pressure due to the re-emitted molecules is 
= 9 = 7 8) / F 
p “pi pRT; (5) = (Vx 2)c,m, ; 
e. « ° ° ec 
Although it is a desirable process from the aerodynamic eiteee 
point of view, specular reflection appears difficult to , 
achieve in practice. The size of interstices between G. = V2RT, 1 
molecules of any ordinary surface is so large that it is m, = mass emitted per unit time tion 
cast 
~ . ° . . - 
For steady-state reflection, m, must equal m,, given by Eq. (4). Hence, for a plane surface, , 
whit 
P, ] C iC, — Uh sin? c Vr C, . ° r / . \ Ras 
- = ST aideakeaiied sin a [1 + erf (UVh sin a)] rs 
(1/2)pU? 2 \U? S ft easi 
or 
Pr l . 7 ba Vr Sil aC, : ; whe 
‘Ay thes rs his {1 + erf (7, sin a)] (10 ve 
('/2)pU? 2M.?7\q 2 Me. & stre 
a gat 
, / ry ‘7 
where ¢,/¢, = /T,/T,. Wi 
for 
> P , , , in < 
Eqs. (8), (9), and (10) summarize the forces to be steady case, the total energy content of streams of rare- de 
. ° - > . 5 
expected per unit area of a plane surface at angle of fied gas is so small that account must be taken of flux nad 
attack a. Ifthe dynamic theory convention isadopted, of radiant energy. Let /, be the net efflux of radiant T 
that an equivalent fraction f of the molecules are dif- energy in unit time per unit area, and let /; and E£, al 
ro . t . . . + Cc t 
fusely reflected, the actual total stresses may be writ- represent similar influx and efflux of energy due to the aay 
ten: gas molecules. oni 
E, = E, — E, 12 
P (9 Ps P; ‘ ; ( tot 
a" S-f2. 3 ts ~ Ai —ws I 
(1/2) pU? (1/s)pU? (1/s)pU? i If flight occurs at night, 
(11) in a 
P re . sd E = BoT 3 
t/('/2)pU? = f[r,/('/2)pU?] ’ AA aot ee" atn 
38 = emissivity = absorptivity of surface 
For nearly all surfaces experimental values of f are ¢ = 5.709 X 10-5 cm. Gm. sec. units, the Boltz- 
greater than 0.90.’ If f is accurately known, Eqs. (11) mann constant 
may be used to obtain more precise values of force .; , : , , ; S: 
ae a P 7 . Daytime flight in the free molecule region will certainly Giv 
coefficients; it is probably a good assumption that f . aan — Ex 
_ Z od . be above any clouds, so that incident solar radiation “q. 
is independent of angle of attack a, so its introduction ees. — : . d 
. , ; 2 must be included in £,. For the upper side of the plate, uc 
would not complicate any integrations of p or 7 over a. solr 
The question of finding the temperature 7, of re- E, = BoT,,' — BE, (sin g cos a’ — sin a’ cos ¢ cos y) e 
5 I : o ? sim 
emitted molecules in diffuse reflection was deferred. E, = radiation density of vertical sunlight 
T, and the surface temperature 7, are fixed by the a’ = angle of attack of surface, positive upward F 
: 8 Bo1 
equality of influx and efflux of energy to and from the from horizontal 
surface. Such equality presupposes a steady-state g = angle of elevation of sun 
process—that all temperatures have had a chance to ¥Y = azimuth angle of horizontal component of ve- wh 
adjust themselves to equilibrium values. For the locity, measured from sun’s azimuth T 
7 wv 
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The energy £; of the incident stream of molecules may be found precisely by an integration such as that used 
The rather complicated result is, for the case of monatomic molecules: 


M.,* sin a\) 
ars ) 


to obtain m; in Eqs. (3) and (4). 
5M. sin a 


E; = : mC,” {(1 + Me) e-uaae “4 4/7 [1 + erf (M., sin a)] ( j 





_. >? sin? a@ 
13° Mo? s ms 


When 1/7. = O, this reduces to 


y~ M,, sina 
f n 


{1 + erf (7, sin a)|t 


E, = mc? = 2m,RT, 


which shows that molecules striking the container have 
kinetic energy */2RT; of the gas itself. 


For \/., greater than about 5, the general expression 
becomes 


E, = m,[(*/2)U? + (5/2)RT]] 


The quantities 7, and -, are found through considera- 
tion of the accommodation coefficient x, defined for this 
case by 


E, — E, = x(E, — Ex) 


where /:,, corresponds to gas at the surface temperature 
T,. For partially diffuse reflection of perfect gas it is 
easily shown that 
x = f[(T% — T)/(T -— Tr)] 

where 7) is the stagnation temperature of the incoming 
stream at 7; A large amount of experimental investi- 
gation has been done on values of x, recently that of 
Wiedmann.’ He finds values between 0.87 and 0.97 
for machined and polished iron and aluminum surfaces 
in air. If no attention is given to possible vibrational 
degrees of freedom of diatomic molecules, it is a fair 
and useful simplification to set x = 1. This gives 
T, = T,,, whence for the O, N atmosphere £, = (3/2)- 
m:RT,, and for the O, Nz atmosphere /, = 2.268 
m,RT,. These last result because the diffusely re- 
emitted gas has zero macroscopic velocity with respect 
to the surface. 

Eq. (12) can now be written in terms of temperatures 
inany particular case. Thus for night flight in an O, N 
atmosphere (J, 2 5): 


BoT ,,* = m, E U? + - RT, — - 


wl Ww 


RT. | (13a) 


t 


Given the conditions of any problem, all quantities in 
Eq. (13a) except 7, are known. Finding 7, = T, re- 
duces to solving a quartic equation. It is most easily 
solved in this case by trial. For a body of any shape a 
similar equation for T,, takes the form (17, 2 5): 


€ 


BoT A, = E U2 + . RT, - % Rre| f mdA, 
and ” - As 
(13b) 
where A, is the total surface area of the body. 7, = 


I’, appears in most equations as the term 


an energy per unit mass exceeding by four-thirds the mean 


c,/¢, = VT,/T, (14) 


The values of aerodynamic coefficients under the 
assumption of diffuse reflection depend usually on alti- 
tude, surface material, and flight 1/7... Under some 
conditions the variation with quantities other than 1/,, 
is small or negligible. Plots of C, and Cp are presented 
for certain bodies; these contain for comparison curves 
of coefficients vs. 7 ., under the assumption of specular 
reflection, where the computation of c,/c; is unnecessary. 
These latter coefficients will be distinguished by super- 
script J, Cp! for example. 

It must be emphasized once more that Eqs. (13) as- 
sumes slightly different forms for daylight flight and 
for the O, Ne atmosphere. Furthermore, below M.. = 
5, account must be taken of the transition 


Ey = ml/g? + (DRT 


to 
E, = 2mRT 


when 1/,, = 0. Reproduction of all pertinent equa- 
tions here is too cumbersome, but, when results are 
listed, the assumptions underlying their computation 
will be given and necessary refinements wi!l have been 
made in that computation. 

The purpose of this section is to give the general 
method for computing aerodynamic coefficients. The 
basic technique consists of integrating the stresses 
given by Eqs. (6) through (10) over the surface of the 
body and dividing by some suitable reference area. 
With diffuse reflection, a further assumption must be 
made on whether the body is composed of a heat con- 
ducting or insulating material. In the latter case, 7, 
will vary with the angle of attack of the surface element, 
according to Eq. (13a); in the former, 7, will be a con- 
stant over the surface, and Eq. (13b) will hold. 

The method is first demonstrated by finding C, and 
Cp for a flat plate at angle of attack a. The shape and 
area of the plate are immaterial, since in free molecule 
flow the action of the stream on each portion is inde- 
pendent of the other portions. The usual convention 
is adopted that C, and Cp are based on the total area 
of the plate; then the forces on unit area consist of 
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appropriate combinations of p(a), r(a@) and p(—a), 
t(—a). The parentheses here refer to the argument 
of the function—the angle at which # or 7 is to be 
evaluated. If the forces Z and D are measured per- 
pendicular and parallel to the free stream, then per unit 
area 


L = [p(a) — p(—a)] cos a — [r(a) + 7(—a)] sina 
D = [p(a) — p(—a)] sin a + [7(a) + 7(—a)] cos a 
If the corresponding expressions are substituted for p 
and 7, 


CoS @ 


CGC = erf (M,, sin a) + 


Me 
Vr sin « cos a (“) (15) 
M. 





9 5 be 
o _ (5) e7 Mot sin? a + 2sinaX 


1 , /n sin? a [c, 
(1 + aad erf (VM. sin a) ae — — ss (2) (16) 


These results are confirmed by Tsien.? For compari- 
son, the same coefficients will now be computed under 
the assumption of specular reflection of molecules. 
The tangential stress 7 vanishes, and 


L' = cos a[2p,(a) — 2p,(—a)] 


C! = 4 sin @ COS @ _ ys sint a 


MoVx 
2 cos @ hc : - 7 
Ga — -+ 4 sin’ a cos «) erf (7, sin a) (17) 
D' = sin a [2p,(a) — 2p,(—a)] 
4sin* a@ _wyiotsinta 
oo a 


@ y/ — 
, M Vr 


2 sin @ ; . 
( 77 _ +4 sin® a) erf ( M. sin a) (18) 


In the specular reflection case, the lift-drag ratio as- 
sumes the particularly simple form 
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Li/D! = CL! Cp! = cota (19) 


These results are confirmed by Snow.‘ Properties of 
the flat plate are plotted in Figs. 1-4. 
It is of historical interest to note that when J/,, be- 


comes extremely large, Eqs. (17) and (18) reduce to 


C,! = 4 sin? a cos a 
Cp! = 4 sin® a 


These expressions for flat-plate lift and drag are men- 
tioned by Zahm.! They are first attributed to Newton, 
who treated the problem of a plate in a uniform stream 
of elastic particles devoid of thermal agitation (equiva- 
lent here to U > c,). They summarize the so-called 
“Sine-squared Law” of flat plate lift, which was once 
employed in attempts to prove the infeasibility of flight. 
When the lift-drag ratio of cot a is considered, even 
these equations do not seem to preclude flight at low 
angles of attack. 

To illustrate analytic integration of surface stress 
distribution to find force coefficients for a body of known 
shape, the drag of a sphere will now be computed. So 


Di = eR? f pi[(7/2) — a] cos a sin a da 
0 


D! t 
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that the example can be as simple as possible, specular 
reflection of molecules is assumed. The results of the 
diffuse reflection case will be given later. The surface 
area of the sphere is dissected into elements of con- 
stant angle of attack as indicated below. 

The circular element of area is formed by rotating arc 
Rda about an axis through the center parallel to stream 
velocity U. If radius R is taken at an angle a to the 
stream, the element has area 27R? sin a da and angle 
of attack [(7/2) — a]. Thisresult is checked by noting 


that 
f 27rR? sin ada = 4rR* 
0 


the total area of the sphere. Assuming specular re- 
flection, the element is subject to pressure 2;[(7/2) — 
a]. Symmetry considerations show that the com- 
ponents of the pressure force normal to the stream can- 
cel one another, giving zero lift. To find the drag force 
the streamwise pressure components 2/;[(7/2) — a] cos 
a must be integrated over the surface [cf. Eq. (G)). 


© 9 ® 
fe - —- e~ Mat cost sin a cos? ada + —— sin a cos ada + 
*y />)p U* xR? MoV 0 Me™ 0 


. ’ [/sin @ cos 
if sin a cos’ ada + 2 f erf (.7., cos a) (R25 *) + 2 sin a cos* a | da 
0 0 1 iw” 


In this case, all integrals can be expressed as ele- 
mentary functions. 

Evaluating and collecting terms and noting that the 
second and third integrals vanish, the final result is 
produced. 


~ I . : : 
Cp! = erf (M.) E r (75) 7 (-) “4 
e- Me? [ 1 2 
oF (Fs) 7 =) | 


This drag coefficient is especially simple for computa- 
tional purposes. It reduces to the value 2.00 for 17, > 
20. . 

As a final example of the generality of the method for 
computing force (and moment) coefficients, it will be 
applied to a case where analytic integration of the pres- 
sure over the surface of the body is impossible. Assum- 
ing specular reflection of molecules, it is desired to com- 
pute Cp! of a typical airfoil at zero angle of attack. A 
symmetrical airfoil of maximum thickness 12 per cent 
chord has been chosen from page 170 of Warner's 
Airplane Design: Performance.* The computation of 
Cp is performed in Table 2 for 4, = 1. The chord 
length is arbitrarily set at 100 units and the spanwise 





length unity. The columns of Table 2 are self-explana- 
tory. Since the airfoil is symmetrical about the chord- 
line and at zero angle of attack, the area AA at mean 
angle a is given by 2 ACh/cos a, as in column 6. The 
pressure on AA in free molecule flow with specular re- 
flection is 2 p; (a). Hence, the streamwise force on 
AA is given by 


AD/(1/2)pU? = [2p,(a)/(1/)pU"] AA sin a 


asincolumn 7. The total of forces in column 7 divided 
by the plan area of 100 square units, becomes the drag 
coefficients Cp’. The value Cp’ = 0.203 may be com- 
pared with about 0.01 in a flow of incompressible, con- 
tinuous air. Greater computational precision may 
readily be had by splitting the airfoil into more chord- 
wise sections and by using more accurate methods of 
tabular integration (Simpson's, Runge-Kutta, etc.). 


Form oF C;, AND Cp FOR SEVERAL SIMPLE BopIes IN 
FREE MOLECULE 


Listed below are some exact or approximate lift and 
drag expressions of six bodies for which they can be 
All forms consist of suitable 


expressed analytically. 
Specular C,! 


combinations of Eqs. (6), (7), and (10). 
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TABLE 2 7 ine il 7 bs 
Numerical Computation of Cp for a Symmetrical Airfoil at Zero Angle of Attack. M. = 1 C 
@ @) @ ® © ®@ 
stati Thick 4A = AD e 
9 a Mean 2p (a) 2ACh (/2)p U2 
nits : ACh a (1/2)p U? COS @ 5) X ® sina 
( 0 = - =e go = — 7 =e 
: 25 3.78 1.25 56.6° 4.67 4.55 17.75 
9'5 518 1.25 29.23 2.69 2.87 3.77 
5.0 7 14 2.5 21.10 2.10 5.36 4.05 
75 8 40 2.5 14.48 1.700 5.16 2.20 
10'0 9 37 2.5 10.98 1.508 5.10 1.460 
15-0 10.70 5.0 6.58 1.334 10.10 1.775 
20 '0 11.49 5.0 1.53 1.190 10.02 0.945 
sO i ae 5.0 5.0 1.048 10.01 0.210 
30.0 12/00 5.0 0.630 1.016 10.00 0.112 in 
10.0 11.60 10.0 —1.145 1.000 20.62 — 0.400 z 
50.0 10.59 10.0 —2.895 1.005 20.04 —1.019 fa 
60.0 9.13 10.0 — 4.185 1.010 20.06 —1.481 
70.0 7 33 10.0 —5.16 1.012 20.08 — 1.830 
80 0 5 OF 10.0 —5.97 1.018 20.11 —2.12 
90.0 289 10.0 —6.70 1.020 20.14 —2.39 
95 0 1 61 5.0 =e 1.027 10.08 —1.312 
100.0 0 24 5.0 = ; 1.030 —_ 10.10 —1.417 ca 
Cp! = 2@/100 = 0.203. 
C, 
and C,! for a particular body in a particular orientation Cp = Cy! + 2WV2/3M)(c,/c:) (21) Vv 
depend only upon the parameter 1/7... C, and Cp M 
based on the assumption of diffuse reflection further To compute 7), and thence c,, the appropriate form of 
require computation of body surface temperature T,,,. Eq. (1b) is easily shown to be 
(A) Sphere.—Based on reference area 7rR?, an 5 
/ r 8BoT,* | ; rT >. 2 V2 - 
- - = * [Zz 
D! 2 " \3pRc “ r " 
Cy! = —— = erf (M..)|2+(——])- 
(1/2)pU*xR* } ‘ , : ae 
. : for night flight in the monatomic O, N atmosphere. 
1 e M | 9 < I pl 
( )] ‘. ; ( ns ) (20) Here the mass flux factor ® can be found by integrat- 
») 1 / 3 ° 7 ~ ° kh 
2M. Va \M. M. ing Eq. (4) over the surface of the sphere. The result 
gives C) 
3 ef 2 . 1 
Co o = 7 a erl (M 4 M.. os = 
Vr 2M . 
» — | | ae 
Computed drag coefficients of the sphere for 17,, up 
to 25 are shown in Fig. 5. Both Cp and Cp! approach 
6 os a a ae ee: 2asM.— ~. 
. , . : C 
| (B) Flat Plate—Eqs. (15) through (19) summarize | 
5 — : ~- a a - the flat plate results. In addition, it can be shown 
| | that 
| | aC. 9 J a 
4 ee — | ae ay - Via [¢, 92 
| = r + (2.9) 
da\.=2.9 MavVr Ma. \co; 
| 
3 22a S | eee re —— — and 
eee ‘ (dC,! da)|. =“ (90 = Pa M.Wr (24) 
Co(/LO km) 
eS a . : = St 
co Representative lift curve slopes appear in Fig. 4. 
(C) Cone with Axis Parallel to the Flow.—Reference 
/ _ < Th 
area is base area, half angle @. 
p 
2 sin 6 ee 
. —M oo? sin? 8 
°° 2 #¢ 6 8&8 10 /2 /4 16 18 20 22 24 Cp! = 7 = € i + 
M., MoV ir 
, ch 
FIG 5 — DRAG COEFFICIENTS OF A - 9 : F ‘ 
A SPHERE VS MM. ( = + 2 sin? 6) [1 + erf (M,., sin @)] (25) ar 
O, N ATMOSPHERE. M.° pl. 
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c FE e . 1 C; |x @/O 
7 sin 0°/ 2M? Ci 
—Mo' sin? @ Ji +5 + - _ -+ (sin ¢ 6) vs (2) x 6 
OM AC; 
[1 + erf (A7., sin @)] (26) 6 
nn ry. ° Co 
Che temperature 7°, is computed from a\ 
— (4 aB0T 3 7 2 - \ 
re ( 2 —— +1)=17,(1 3 M ? (27) S ge 
5pRc& 3 ee 2 
° . _ ” Co 
in the monatomic O, N atmosphere. The mass flux 
factor ® can be found directly from m,, 
—Moo? sin? 0 mae : : CO 2 ¢ 6 8 10 12/4 16 18 20 22a 
ge=e"* "+ MoV rr sin 6 [1 + erf (.V/,, sin 6)] Mu 
: eae . , FIG 6 - DRAG COEFFICIENTS OF 
(D) Double Wedge Airfoil.—For a doubly symmetri- CIRCULAR CYL/INDER NORMAL 
aire ae ee ee ae ae TO STREAM. / N O 
cal airfoil of leading edge angle 26 at zero incidence, pi Avan ~ aut = pe ‘ R-# 
9 
c ee (; es ) e M ~? sin? 6 , 
= = 
MoV x cos 8 C1 4 sin @ tan ‘) —Maotsin?d 
= — é 7 t 
D Mav . 


om , 
va (2) sin @ tan 6 + 2 tan 0 (1 + a) x 


ae + tat Oemelet i. ae OD 
erf (M. sin 6) (28) sin an er! (/,. sin (Zt 


M..? 


Lift and drag coefficients at angle of attack a and lift curve slopes can be found by proper combination of flat 
plate results. For example, in the doubly symmetrical case of leading edge angle 26 


l 6) — 6 
CGC = —. on — erf [.7,, sin (a + @)] + rae be erf [.\7., sin (a — @)| + 
Vr . 
ve (6 ) sin (a + @) cos (a + 4) +> —— = (2) sin (a — 0) cos (a — a) (30) 
co Ci 


(£) Right Circular Cylinder with Axis Perpendicular to the Stream.—Reference area is plan area of cylinder. 


_{2n + 3 2n + 1 
M. ot (? _) ) Pas r( t 3 ) 


9 * 
Cp! = — —1)" hen Seer a —1) re" 
M a ! I'(n + 2) M > ( n! I'(n + 2) * 


on= 0 nN: on= 


2n+ 3 _{2n + 1 
in gp emurhe®), HH) 


M (—1) 
5) 4, n! | C(n aa 3) '(n fe 2) 


rp (22+ ') 
Ps M = 2 


(—1) 32) 
a n! C(n + 2) ° 


r (> + ') 
M.”" 2 rV/ rC, I 
‘ + + (ar. + ) 


= : om 59 
Mo , | n! T(n + 1) 4M .C; 


_{2n + 1 ,{2n + 1 
pet ( 2 ') v2! ( 2 ) 


tBol 4 l 3/5 Pn ) : | 
“ =T7,|-M,? — = )] (—1)” + M.? —1)" 
pRC,; 3 T 4 (; T RF | n! T'(n + 1) a> \ n! C(n + 2) \ 


Cp 


Surface temperature and C,/C; can be computed from 


(F) Circular—Arc Ogive.—This is one-half the figure produced by rotating a minor arc of a circle about its 
chord. If the reference area is the base area, the body is described by its nose angle 26. The expressions below 
are based on the assumption of moderately small #6. C,/C, can be approximated by its value for one side of a flat 
plate at a = 6/2. 
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where V = (1 + cos 6)/2. 


NUMERICAL RESULTS 


(A) The flat plate and ogive drag coefficients can be 
used to estimate drags of bullet-shaped bodies. Con- 
sider a bullet consisting of ogive nose on a circular 
cylinder of length / and radius R. As a numerical 
example, let R = '/2 cm., / = 3 cm., and nose angle 
20 = 20°. If MW. = 2 at 120 km. altitude, the flight 
speed U = 1,021 m. per sec. and Cp = 8.20. 

In the hypothetical specular reflection case Cp’ = 
0.229. Increase in \/,, greatly reduces the drag co- 
Thus, under the same conditions at 120 km. 
= 10, U = 5,100 m. per sec., 


Cy = 3.77: Cp! = 0.0650 


efficient. 
with 1/7, 


These drag coefficients in the free molecule region are 
to be compared with roughly 0.20 at sea level and Mach 
Number 2 and 0.05 at Mach Number 10. To contrast 
actual drag forces, they would be 4.34 X 10° dynes 
at sea level, 12.29 dynes at 120 km., and 1.181 dynes at 
180 kin. 
less in the free molecule region, unless flight speed is 
much greater. At 120 km. and jA/.. = 10, the drag has 
risen to 140 dynes, and it would be several thousand 
dynes above 1/.. = 20. The drag could be made truly 
negligible by using a specularly reflecting material, in 
which case it would drop to 0.343 dynes at 120 km. and 
MM, = 2and 2.54 dynes when IJ, = 10. 

(B) Rather uninhibited speculation about the pos- 
sibility of aerodynamically sustained level flight at free 
molecule altitudes leads to the following conclusions: 
The best lifting surface is the nearest thing structurally 
possible to a flat plate. Adequate L/D can only be ob- 


Thus, the force on such a small object is vastly 
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tained from a specularly reflecting material. The wing 
loading would have to be almost ridiculously small. 
As an example of the necessary orders of magnitude, a 
500-kg. airplane with 10,000 sq.m. (2.47 acres!) wing 
area will fly at 120 km. as slow as JJ, = 10. The lift 
coefficient would here be 0.104 and the drag force a 
mere 69 kg. If flight were westward above the equa- 
tor, the upward centrifugal force would only be about 


190 kg. Moreover, this airplane would also go to 
higher altitudes in the free molecule region; for ex- 


ample at 180 km., the minimum 1/7, for practical flight 
is roughly 18, at C,’ = 0.4. None of these aircraft 
could fly in the extremely rare gas above 220 km. 
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Some Aeroelastic Properties of Swept Wings 


S. I. PAI* anp W. R. SEARS? 


Cornell University 


ABSTRACT 


The problem of the lift distribution on an elastic wing is re- 
examined, especially from the standpoint of the swept wing. It is 
shown how, using matrix notation, the problem can be set up 
and solved in a rather general class of situations. 

As examples of this procedure, certain typical numerical cases 
involving a sweptback and a sweptforward wing are treated. 
The effects of sweep on certain aeroelastic phenomena, such as 
torsional divergence and aileron reversal, are confirmed. The 
case of dive pull-out is a typical one involving several variables 
related through the dynamics of the airplane; the aeroelastic 
effects may be large. 

Finally, a case involving an airplane longitudinal oscillation, 
together with wing bending, is investigated. Here it becomes 
practically necessary to assume the mode of wing bending de- 
flection. The results for this numerical case indicate that elas- 
ticity modifies the short-period longitudinal oscillation in such a 
way that the effective static stability and damping in pitch are 


reduced. Unstable oscillations may occur in extreme condi- 


tions. 
SumMaryY oF NortaTion{ 
a = local slope of lift curve for the profile at s, cor- 
rected for Mach Number 
(t = operator giving incidence distribution required to 


produce the lift-coefficient distribution c;(s) 

h/2 = value of y corresponding to s = s,; hence, the pro- 
jection of s,; on the transverse axis; for complete 
wing, the half-span 

local chord measured parallel to flight direction 

r mean chord for definition of Cy 

local coefficient of lift, based on g and ¢ 

local coefficient of moment about aerodynamic 


m 
center, based on g and c? 

te = local normal-load coefficient based on g and c; 
taken equal toc 

Ci = rolling-moment coefficient 

Cr = lift coefficient of airfoil, based on g and S 

Cu = moment coefficient of airfoil, referred to airplane 
c.g., based on g, S, and z 

d = distance of local center of mass forward of elastic 


axis, measured perpendicular to elastic axis 
e = distance of aerodynamic center forward of elastic 
axis, measured perpendicular to elastic axis 


EI = bending stiffness distribution 

[ E] = ({H][ce] — [F][c])[S] 

f(s,o) = [in Eq. (2.11)] any function of sand o 

fpr = S (sp, Cr), (p, ae Rs 6 sie n) 

[f] = the square matrix whose members are the values 


Sor 


F (s,o) = tF(s, 0) + yF,(s, a) 


Presented at the Annual Summer Meeting, I.A.S., Los Angeles, 
July 14-16, 1948. 

* Visiting Professor, on leave from National Central Univer- 
sity, Nanking, China. 

t Director, Graduate School of Aeronautical Engineering. 

t Notation used only in the Appendix has been omitted here 
for brevity. 
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Fi(s, a) 
F,(s, a) 


g 
g(a) 


GJ 
H(s, o) 
Ai,(s, o) 


H,(s, a) 


“ys 


[S] 


ea™ 
+ 


x1/, 
xI/, 


= 


influence function giving zs, at s due to unit normal 
load at o 

influence function giving 8 at s due to unit normal 
load at ¢ 

acceleration due to gravity 

[in Eq. (2.11)] any function of ¢ 

g(or), (ry = 0,1, 2,... 2m) 

the column matrix whose members are the values 
&r 

the diagonal matrix whose members are the values 
gr. (ry = 0, 1,2,...2) 

torsional stiffness distribution 

tH,(s, ¢) + yH2(s, o) 

influence function giving zs at s due to unit torque 
applied at ¢ 

influence function giving 8 at s due to unit torque 
applied at o 

moment of inertia of airplane, in pitch, about its 
c.g. 

moment of inertia of airplane, in roll 

unit matrix, whose diagonal members are ones and 
all other members zero 

((F] — [A}[d})[ 5S] 

mass distribution: 

mass of airplane 

stability derivatives giving pitching moment per 
unit w, 6, for rigid airplane 

number of subintervals into which the range of 
integration 0 < o < s; (or 0 < 9 < 8/2) is di- 

hence, the order of the matrices used is 


mass of wing per unit s or ¢ 


vided; 
n+ 1 

local load factor due to local normal acceleration 
and component of gravity 

linear operator giving cc; corresponding to any 
given incidence function 

dynamic pressure of flight 

radius of curvature of flight path 

distance measured along elastic axis from one end 
of airfoil 

value of s at other end of airfoil 

the values of s at the end of the rth subinterval— 
ie.,7 A 

wing area 

integrating constant for the 
For the trapezoidal rule, the values of Sy, are 
A/2, A, A, ... A/2; for Simpson's rule 4/3, 
4A/3,2A/3,...44/3, 4/3. Ais the length of a 
subinterval 

the diagonal matrix whose members are the values 
S, (ry = 0,1,2,...#) 

time 

velocity of flight 

amplitude of wing bending deflection at 5, 

velocity of airplane c.g. in z direction (body axes) 

airplane gross weight 

distance forward from y axis (through airplane 
c.g.) to elastic axis 

distance forward from y axis to aerodynamic center 

distance forward from y axis to three-quarter- 
chord point of profile; * x + ey — (c/2) 


rth subinterval. 
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= distance forward from y axis to center of gravity 
of profile; ~ x + yd 
y = projection of s on transverse axis 
elastic deflection of elastic axis normal to the plane 
of the airfoil (positive downward) 


Zw, Zg = stability derivatives giving vertical force per unit 
w, 6, for rigid airplane 

a(s) = local absolute angle of attack 

a) = local angle of attack due to elastic deflection 

a’ = local absolute angle of attack of rigid wing, in- 


cluding control deflections and airplane linear 
and angular velvcities 
3 = elastic deflection of airfoil in twist about the elastic 
axis (positive nose-up) 
cosine of local angle of sweepback of elastic axis 


= aileron deflection 


a 
” = same as y 
curvature of the elastic axis, in the plane of the 
wing 
r = complex exponent 
. = Xr 
u = JI/pSb 
v = ratio of wing weight to total airplane weight 
= 25,m/SW 
= sine of local angle of sweepback of elastic axis 
) = mass density of air 
a = same as Ss 
T = SIU/pSU 
o(s/s1) mode of fundamental, wing bending at zero air 
speed. Note: ¢(1) =1 
v = frequency of fundamental wing bending at zero 
air speed 
y’ = @T 


(I) INTRODUCTION 


A NUMBER OF AEROELASTIC phenomena of airplanes 

have been recognized for a long time—the 
classical cases of wing torsional divergence and aileron 
reversal are most familiar. In addition, the more gen- 
eral nature of the aeroelastic problem has been pointed 
out by several authors.'~* 

The present paper deals with the particular case of 
an elongated wing having sweepback or sweepforward. 
It appears that aeroelastic effects are exaggerated by 
sweep, especially because bending of a swept wing di- 
rectly affects the lift distribution. * 

The paper is divided into three main parts. 
first, consideration is given to the fundamental aero 
elastic problem, the calculation of the aerodynamic 
load distribution on an elastic airfoil. In the second, 
the methods developed are applied to certain typical 
problems, the results of which exhibit the characteristic 
Finally, a typi- 


In the 


aeroelastic behavior of swept airfoils. 
cal dynamic problem involving both a rigid airplane 
and its elastic wing is attacked. 
(II) THe AERODYNAMICS OF AN ELASTIC AIRFOIL 
The considerations of this section are restricted to 
cases of “elongated’’ aerodynamic surfaces—i.e., those 
* The authors know of two earlier, unpublished, investigations 
of some aeroelastic phenomena of swept wings; one by J. W 


Miles, at Northrop Aircraft, Inc., the other by F. W. Diederich, 


at Langley Field 


that have one dimension, nominally spanwise, much 
greater than the chords, and for which the basic as- 
sumptions of the Prandtl wing theory are substantially 
valid. For such wings one can assume that the aero- 
dynamic load distribution is completely determined by 
the ultimate distribution of incidence along the span 
and, conversely, that the ultimate elastic deflection at 
any station is uniquely determined by the spanwise 
distribution of normal force and torque. 

Let s be the distance measured along the airfoil in 
the direction of its principal elongation. In view of the 
assumptions already made, it seems permissible to 
assume that the airfoil possesses an elastic axis; then 
s w Il be measured along this axis (Fig. 1). The inci- 
dence a(s), at any condition of flight, can be expressed 


as 


als) = @*? + eg’ (2.1) 


where a’ includes angle of attack, built-in twist, and 
effective incidence due to control deflections, angular 
and linear velocities, etc., for the rigid airfoil; and where 
a’ denotes the elastic effect. There seems to be no 
essential value in trying to write out a general expres- 
sion for a’; the variety of terms it may include seems 
almost endless, since the shape, location, and function 
of the airfoil under consideration have not been speci- 
fied. Specifically, however, a(s) is supposed to be 
measured in the flight direction and does not include 
induction effects of the airfoil itself. 

The elastic incidence a‘? depends, in a swept wing, 
on both the twist of the airfoil about its elastic axis and 
the bending of that axis; the equations that control 
these phenomena are assumed to be the curved-beam 


formulas:* 6 


[/:1(z;; — Bx)],, = normal loading ats (2.2) 
[GJ(8, + 2,x)], = torque loading at s (2.3 


where the bending deflection z and the normal loading 
are taken positive in the same direction, say downward, 
and similarly for the twist 8 about the elastic axis and 
the torque loading. The subscripts denote differentia- 
tion. Now the solutions of Eqs. (2.2) and (2.3) can 
always be expressed in terms of four influence fune- 
tions, Fi(s, 7), Fo(s, o), Ai(s, o), and H2(s, o), which 
give, respectively, the values of z, and 6 at s due to unit 
normal load and unit applied torque ato. In some cases 
it may be possible to determine these influence func- 
tions experimentally. Furthermore, the normal load 
and torque are supposed to be the results of aerody- 
namic loads, inertia loads, and weight, so that, symbolic 
ally, the important quantities 2, and 6 can be expressed 


as follows: 


F,(s, ){ —qe,c + mNg}do + 


0 


| M(s, a) {glcnce + Cnc?y] — mNegd}do (2.4) 


Z,(5) = 
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B(s) = Fé F2(s, ¢){ —gcenc + mNg}do + 
0 
f Fy(s, 7) {gence + Cnc?y] — mNegd}do (2.5) 
0 


Here c, denotes the local normal-force coefficient, 
which will be taken equal to the lift coefficient c;. The 
terms mNg and mNgd represent the weight and inertia 
loadings, where JN is the local normal “‘load factor’ at 
o—i.e., the absolute acceleration of the airfoil in the 
—z direction increased by the appropriate component 
of gravity. 

The elastic incidence a“? is given by 


a = y8 + és, (2.6) 

Now, to determine the load distribution and the de- 
flection of the elastic airfoil, Eqs. (2.1), (2.4), (2.5), and 
(2.6) have to be put together and c, and c,, determined to 
correspond with a(s). But since c,, has been referred 
to the aerodynamic center of each profile, it is independ- 
ent of a and can be omitted from this statement. 
Symbolically, introducing the notation @[c,(s)] to de- 
note the incidence distribution a(s) required to produce 
the lift distribution c,(s), one can now write 


a(s) = @[e,(s)] 


a+ fi H(s, 0) {qlewe + ¢pc?y] —mNegd}do + 
0 


f F(s, a){—qec + mNg}do (2.7) 
0 


Eq. (2.7) is, of course, an integral equation for the 
unknown function c,(s). The functions a”, 6, €, Cms 
y, m, N, and d are supposed known, as well as the elas- 
tic influence functions 7 and F. Some of these depend 
only on the geometry of the airfoil, its mass distribu- 
tion, and the location of its elastic axis, while a, 
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Fic. 1. Diagram showing notation 


Cm, and N depend also on the control deflections and 
airplane motions involved. Actually, these are not 
independent but are coupled through the dynamical 
equations of the airplane—or some simplifying rela- 
tionship that replaces them—and these generally in- 
volve integrations of c,(s). It is not particularly 
helpful to point out that this problem is a complicated 
one, but it is desirable to recognize these interdepend- 
encies so as to avoid errors. 

The practical procedure to solve the integral equa- 
tion, Eq. (2.7), is to adopt matrix notation and to ac- 
cept an approximate numerical answer. Before pro- 
ceeding in this way, it may be advisable to discuss the 
operator @[c;]. 


The Operator @ 
The simplest choice of this operator is that corre- 


sponding to strip theory: 
Q@.= i/a 2.8) 


where a is the slope of the lift curve, presumably taken 
to be a constant for the entire span. Usually the two- 
dimensional value is not used for a but rather a reduced 
value appropriate to the aspect ratio of the wing. Ob- 
viously, this approximation is a crude one, but there are 
many investigations for which more refinement is un- 
necessary. 

To account properly for induction, at least for wings 
of small sweep, the familiar Prandtl wing theory should 
In this case the operator @[c,] stands for 


th a Oe 4 L [ea dy 
a(s) Str dn y — 7 


Here y and » represent the projections of s and o on a 
transverse axis, and the integration must be carried 


be used. 


(2.9) 


over the entire span. 

For highly swept wings, certain refinements or ex- 
tensions of the Prandtl theory are available.6-* Of 
these, the ‘‘L-method”’ of Weissinger’: ® seems to be the 
best for practical use, even though its justification is 
surely more empirical than theoretical. In any case 
it will be assumed in this paper that the operator is 
linear, as is the case with all of the approximations men- 


tioned here—that is, 


ale, + c,] = Alc, ] + alc, ] 2.10) 


In a later paragraph it will be shown how all of these 
operators @ can be expressed in matrix notation for 


practical work. 


Formulation Using Matrices 

If the various integrations that appear in the equa- 
tions above are approximated by finite sums, using a 
convenient integrating formula such as the trapezoidal 
rule or Simpson’s rule, the equations can be written in 
matrix form. The basic approximation and defini- 
tion of this process are the following: 
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MMs) = f "Ki, Mtoe « DIG, eeS, (210 
0 r=0 


or 
h(sy) =h, ~ e602 (p = 0,1, 2,...m) (2.12) 
r=0 
or 
th} = [fJ[S]tg} (2.13) 


where |/} and {g} are column matrices whose members 
are the values h, and g,; [f] is the square matrix whose 
members are the values f,,; and [S] is the diagonal 
matrix whose members are the values S,. 


Thus one proceeds to define a column matrix} } anda 
diagonal matrix [] corresponding to each function of 
s, and a square matrix [ ] corresponding to each func- 
tion of s, ¢. Similarly, the linear operator @, for any 
of the approximations suggested above, is replaced by a 
square matrix [ @], defined by 


fa} = [@]{c,} (2.14) 


For example, the matrices [@] corresponding to strip 
theory, Prandtl theory, and Weissinger’s approxima- 
tion are given in the Appendix to this paper. 

The fundamental aeroelastic equation, Eq. (2.7), can 
now be stated in matrix form: 


[@}fe} = {a} + LF] —glelSI{ea} + 


gIN][S]tm}) + [A] (glee][Slte.} + gle*}[v] Xx 


[S] tm} — gld][N][S]{m}) (2.15) 
or 
[@]{ey = ta} + gM] [ec] [y][S] len} + 
g([K][N]{m} + @g[E]{e.} (2.16) 


where [Z] denotes ([H][ce] — [F][c])[S] and [AK] de- 
notes ([F] — [H][d])[.S]. 


The solution of Eq. (2.16) can be indicated by a con- 
venient symbolism : 


tes} = ([@] — g[E])-“(fal’} + 
Q(A} [e*][y][S} tem} + g[K][N] {m}) 


(2.30) 


where ({@] — g[£]) ' denotes the reciprocal matrix to 
[@] — g[E]—ie., their product is [J]. The problem 
of finding the reciprocal matrix is discussed in several 
textbooks.? 
q, the dynamic pressure, is a given constant, and there 
are more complicated procedures that yield additional 
information when q is a variable parameter. Since 
[@] depends upon the Mach Number, these more 
elegant procedures may not be particularly useful in 
cases relating to high-speed flight. 


There are simple methods for cases where 
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Fic. 2. Wings used for numerical examples. 





Total Lift, Etc. 


The lift-coefficient distribution {c,} thus calculated 
is now available for integration to determine the, total 
lift on the airfoil, the aerodynamic moment about any 
desired axis, and other useful data. 

It may be advisable to review the meanings of the 
terms in Eq. (2.17). 

(1) The following matrices are given by the geom- 
etry and elastic properties of the airfoil and its sup- 
porting structure: [/], [//], and [A]. 

(2) The following depend on the geometry alone: 
[c?], [vy], and [NV]; [NV] expresses the normal accelera- 
tion at s in terms of the acceleration components of the 
airplane (at its center of gravity) and its angular ve- 
locities and accelerations. 

(3) {m} represents the mass distribution along the 
airfoil. 

(4) {cm} and {a} depend on the aerodynamic con- 
figuration of the rigid airfoil—i.e., the size and effective- 
ness of its control organs, their deflections, and the air 
plane’s six components of linear and angular velocities. 
Gust velocities would also affect these matrices. 

Obviously the integrated quantities, such as Cj, cal- 
culated from {c,} in Eq. (2.17), will have analogous 
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TABLE 1 
Properties of Example Wings 





Area 271.36 sq.ft. 
Chord (measured parallel to axis of symme- 

try) 8.48 ft. 
Chord (measured perpendicular to leading 

edge) 6.00 ft. 


Elastic axis located at 33 per cent chord from 
leading edge 

Center of gravity located at 43 per cent chord 
from leading edge 

Bending stiffness, EI 

Torsional stiffness, GJ 


23.7 X 108 Ib.ft.? 
2.39 X 10° lb.ft.? 


Span 32.0 ft. 

Wing weight To be specified 
Profile Symmetrical 
Built-in twist Zero 
Sweepback: Wing No. 1 45° 
Sweepforward: Wing No. 2 45° 





terms depending upon the same variables. In par- 
ticular, for a given airplane at a given speed and alti- 
tude, these will depend on the control deflections and 
the linear and angular velocities and accelerations. 
Strictly, these are all related through the airplane’s 
equations of motion, and the forces appearing in the 
equations of motion are the integrated quantities that 
are being calculated. It is clear that the situation is 
complicated, and that some sort of simplified attack 
on individual problems must be made. 


(III) Examples: STEADY FLIGHT 


To illustrate the application of the results above, 
certain typical calculations will be carried out here. 
As simplified but rather extreme examples, attention 
will be focused on two hypothetical airplane wings, 
shown in Fig. 2. These wings have the physical char- 
acteristics* as shown in Table 1. 


Lift Distribution 

The lift-coefficient distribution in steady unacceler- 
ated flight at any angle of attack a is calculated di- 
rectly from Eq. (2.17), putting {a} equal to a{J}, and 
im and [LN] equal to zero: 

{c.} = ([@] — g[E]) {Tha (3.1) 

For Wings No. 1 and No. 2, because of their uniform 
dimensional and elastic properties, it is relatively easy 
to determine [EZ]; its members are the values E,,S, = 
E(s,, @)Sp, where 


Me gh a 2 we. etwe 
“Ss, o) = a oe = O's o AY 
GJ 2EI 
: (3.2) 
Ce fc | 
o =e -(2s~s), Cf sXe 
GJ 2EI 


For the aerodynamic matrix [@], the form corre- 
sponding to Weissinger’s approximation is chosen, tak- 
ing the value 27 for the profile constant (including 
Mach-Number effect). The results of this calculation 


* These are, essentially, sweptback and sweptforward ver- 
sions of the wing considered in reference 10. 


for Wings No. 1 and No. 2 are shown in Fig. 3. It was 
necessary to use widely different values of g, the dy- 
namic pressure, for the sweptforward and sweptback 
wings, in order to obtain reasonably comparable re- 
sults, since the aeroelastic effects are so much more 
violent in the former, as will be mentioned again below. 

It is seen that the sweptback wing shifts its load 
toward its center at high speeds, while the sweptforward 
wing shifts its load outboard. These effects, which are 
probably generally known, result principally -from the 
bending deflection. In the sweptforward case the 
bending is in the direction to reinforce the usual tor- 
sional-deflection effect, while in the sweptback case 
they counteract one another. 


Torsional Divergence 


The speed at which the air loads exceed the elastic 
restoring tendencies of the wing structure is given by 
the value of g, for which c; increases without limit for 
any given a, in Eq. (3.1). ‘This can only happen when 
the terms of ({@] — g[£])~ become infinite, and the 
condition for this is that the determinant of [@] — 
q[E] vanishes—i.e., 


[@] — g[E]| =0 (3.3) 


This condition yields a polynomial for g, whose de- 
gree is m + 1, the order of the matrices. The diver- 
gence speed has been calculated for Wings No. 1 and 
No. 2, with the following results: 

Wing No.1: (no positive root). 

Wing No. 2: g = 625 Ibs. per sq.ft. (indicated 
speed = 495 m.p.h.). 
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Fic. 3. Effect of elasticity on lift distribution. Solid curves: 
rigid wing. Dash curves: (upper) Wing No. 1 at g = 640 lbs. per 
sq.ft.; (lower) Wing No. 2 at g = 100 Ibs. per sqft. 
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Fic. 4. Aileron effectiveness and reversal. (Note that the 
calculation is valid only up to a certain subsonic Mach Number.) 


This confirms a known result—namely, that the di- 
vergence speeds of sweptforward wings are often ex- 
tremely low, while the phenomenon of divergence can 
be completely eliminated by the use of sufficient sweep- 
back. 


Aileron Effectiveness and Reversal 


The lift distribution {c,} due to aileron deflection 
5, is obtained from Eq. (2.17) when {a‘”’} and f{c,,} 
are the distributions appropriate to aileron deflection, 
both proportional to 6,, and [N] is the local load fac- 
tor due to rolling acceleration. The members of the 
diagonal matrix [NV] are the values AN(s,), p = 0, 1, 2, 


...n, where 
AN(s) = (qgSb/gI,,)vC, (3.4) 


The rolling-moment coefficient C, at zero rolling ve- 
locity is obtained by appropriate integration of c,; 
thus it is given by a formula containing a term propor- 
tional to 6, and another proportional to C,, and, finally, 
in the form (dC,/d6,)64, where dC,/d6, depends on both 
qand /,,. The vanishing of dC,/dé, defines the aileron 
reversal speed; this is determined by another poly- 
nomial equation for g, independent of J,,. 

The result of numerical calculations for Wing No. 1, 
with ailerons as indicated in Fig. 2, is presented in 
Fig. 4. The reversal speed is well outside the range of 
Mach Numbers for which the calculation is valid. A 
similar calculation for Wing No. 2 results in a value 
of 2,340 m.p.h. for the indicated reversal speed. Na- 
turally, this value is without practical significance, even 
without regard to the Mach Number, since this wing 
encounters torsional divergence at a much lower speed. 
In fact, the same elastic behavior that produces the low 
divergence speed results in a high reversal speed. Here 
the bending under load assists the aileron, in opposi- 
tion to the ‘‘trim-tab”’ effect of the aileron torque. 

It is interesting to note that the wing bending de- 
flection due to rolling acceleration—the term arising 
from [N]—assists the ailerons of the sweptback wing 
as shown in Fig. 4. For a sweptforward wing this con- 
tribution opposes the aileron. The calculation for 
I,; = © also represents the case of a model wing in a 
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wind tunnel, and it is clear that such wind-tunnel re. 
sults require careful interpretation, especially if the 
model has appreciable elasticity. 


Dive Pull-Out 


The tendency of a sweptback wing to unload its tips 
by elastic deflection under load—the same effect that 
eliminates the phenomenon of wing divergence—is 
clearly a destabilizing one in airplane dynamics. This 
inboard shift of the load is also a forward shift and may 
result in a dangerously far forward position of the ef- 
fective aerodynamic center. This is a case typical of 
the profound effects of aeroelasticity on the airplane's 
stability characteristics. 

On the other hand, the increment of angle of attack 
that produces this effect is presumably accompanied 
by an increment of normal load factor, which introduces 
a counteracting aeroelastic effect. Once more this is a 
typical case, for it illustrates the necessity of consider- 
ing all the variables in Eq. (2.17) and emphasizes the 
fact that these variables are actually related through 
the dynamics of the airplane. 

The problem must be more definitely specified. Con 
sider the case of an airplane in a steady pull-out at 
practically constant speed, U. The air loads on the 
elastic wing will result from angle of attack, angular 
velocity in pitch, and load factor—i.e., centrifugal 
These variables are all related to one another, 


force. 
for, if R denotes the radius of curvature of the flight 
path, 

AN = U?/gR = qSAC,;/W (3.4 
and 

6 = U/R = (qgS/UW) AC, (3.5 


It is assumed that the tail provides the balancing 
pitching moment to maintain this flight condition 
during the pull-up. The increment of pitching-moment 
coefficient due to wing, AC,,, for a given value of A.V 
or, according to Eq. (3.4), the increment ACy, for a 
given increment AC,-—measures the stability of the air- 
plane in the pull-out. Instability is characterized by 
positive values of this ratio AC,,/ \C;, while extremely 
large negative values would indicate a loss of maneuver 
ability. 

Again Eq. (2.17) is the basic relation used. In this 
problem c, is again zero (approximately), while a’ 
Aa and the distribution 
The inertia term 


includes the constant value 
—x./6/U due to pitching velocity. 
has the constant value of AN, in terms of ACy as in 
Eq. (3.4). Thus, after integration, one has expressions 
for AC; and ACy, as linear functions of Aa and AC;: 


AC, = ki Aa + ke AC, 
ACy = ksAa + Rk, AC, 


The result, in this notation, is, of course, 





Opp 
vah 
§ = 
plat 
vali 
wit 
effe 
rem 
effe 

F 
poit 
Thi 
ing 


voly 


nel re- 
if the 


its tips 
‘t that 
1ce—is 

This 
d may 
the ef- 
ical of 
ane’s 


attack 
panied 
duces 
lis iS a 
isider- 
es the 
rough 


Con- 
ut at 
ym the 
igular 
‘ifugal 
other, 

flight 


cing 
dition 
ment 
AN 

for a 
e ailt- 
dl by 
-mels 


uver 


1 this 

a : 
ution 
term 
as in 
sions 


my 


(3.6 





AEROELASTIC PROPERTIES OF 


TABLE 2 
Dive Pull-Out—Wing No. 1 


Wing Weight 


_ Wing Loading q | 
Lbs. per sq.ft. | 
I 15 640 4 
2 | 15 640 8 as 
3 15 640 , 8 outboard 
| 4 inboard 
{ 15.6 640 8 
5 45.6 640 45.6 
ACy/ AC, = [(1 — kz) /Rilks + Ra (3.4) 


All of these &’s are calculated by straightforward 
matrix multiplication involving the fundamental ma- 
trix ([@] — g[E])-". 

This calculation has been made for Wing No. | for a 
number of cases. In addition to the data given for this 
wing in Table 1, the location of the airplane center of 
gravity has been taken at the aerodynamic center of 
the wing; thus the wing would contribute nothing to 
the airplane static stability in level flight if the wing were 
rigid. 

The results are collected in Table 2, where values of 
AC,/ Aa and AC,,/ AC; are tabulated for several differ- 
ent wing loadings and wing weights. For comparison, 
the same ratios are given for a geometrically identical 
rigid wing. It is seen that aeroelastic effects are al- 
ways destabilizing, especially when the wings are rela- 
tively light. In Case 3 the wing weight is twice as 
great in the outboard half of the wing panel as in the 
inner. This configuration, which might represent a 
wing carrying heavy stores near its tip, produces the 
most stable airplane in the pull-out. Case 5 repre- 
sents an idealized flying-wing airplane, since all of its 
weight appears in its wing. 

It is important to point out that the two aeroelastic 
effects, due to air load and inertia load, respectively, are 
opposite in sign. For all five cases in Table 2, the 
value of OC,/OC; calculated for the elastic wing with 
§= 0 = AN (as if the wing were attached to an air- 
plane of infinite mass) is 0.0235; the corresponding 
value of O0C,/Oa@ is 3.254. Comparison of these values 
with those tabulated in Table 2 shows the combined 
effect of inertia and angular velocity. Again we are 
reminded that wind-tunnel results do not include these 
effects and may require careful interpretation. 

From these calculations it is seen that the ‘““maneuver 
point’’ is shifted forward appreciably by elastic effects. 
This destabilizing tendency may be minimized by add- 
ing weight to the wing, especially near its tips. 


(IV) EXAMPLE INVOLVING OSCILLATIONS 


Che last example of Section III obviously is one in- 
volving accelerated, rather than ‘‘steady”’ flight, but 


Plan Form as Wing No. 1 
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Rigid Wing of Same _ 
Wing No. 1 





AC / AC; | 
| 
| 
| 





ACL / Aa ACs / Ae ACu/ACz 
4.040 —0.0238 3.400 0.0183 

~ 4.040 ~ 0.0238 3.590 0.0168 
4.040. —0.0238 3.520 0 0129 7 

~ 4.010 —0.0078 3.360 0.0211 | 
4.010 —().0078 3.870 0.0168 


it has been included in that Section because the ac- 
celeration involved is a steady one. Consideration 
will now be given to an aeroelastic phenomenon that 
Actually, this 


is not novel, since every case of airfoil flutter is an aero- 


involves an oscillatory type of motion. 
elastic oscillation. In this case, however, it is pre- 
suined that rather slow vibrations may occur, involving 
only one of the wing vibrational modes, say bending, 
and in addition the airplane’s longitudinal oscillatory 
motion—the ‘‘short-period oscillation.”’ 

One's interest may be directed to this case in view of 
the favorable effects of relatively large wing weights in 
Section III. Is there not a certain risk in increasing 
the wing weight, for one is thereby reducing the na- 
tural frequencies of the wing? 


Equations of Motion 


Assuming that the oscillations will be slow, one may 
neglect nonstationary aerodynamic effects and proceed 
to apply the equations of the preceding Sections with 
minor alterations. The airplane velocity U will be 
taken to be constant; the motion of the rigid fuselage 
and tail are described by the normal velocity w of the 
airplane c.g. and the pitching velocity 6 of the fuselage 
and tail about this c.g., as in airplane dynamics; these 
are referred to body axes. 

The differential equations of wing deflection, Eqs. 
(2.2) and (2.3), now become partial differential equa- 
tions, involving both s and the time, ¢. The problem, 
while certainly possible of solution by numerical ap- 
proximations, is now too complicated to be attractive. 
Nevertheless, there is a powerful tool of the structural 
engineer that seems particularly appropriate for use 
here—namely, the Rayleigh-Ritz approximation.* This 
method involves the assumption of the modes of elastic 
According to the 


deflection or a series of modes. 


theory, an approximation to the fundamental fre- 


* AH. Flax points out that the Rayleigh-Ritz theory is not 
strictly admissible here, since self-adjoint differential equa 
tions are assumed in that theory. The partial differential equa 
tions of the wing with sweep, including air loads, are not self 
adjoint, so that an extension of the theory seems to be required 
Intuitively, at least, one does expect.the method to yield good 
approximations even in this case. If the modes assumed are 
good approximations, the results will surely be dependable. 
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quency of vibration is obtained even if the choice of 
mode is poor. 

One might legitimately ask why it is permissible to 
assume a mode of wing deflection at this point, whereas 
the deflection curve has always been calculated in ma- 
trix form heretofore in this paper. The answer is that, 
whereas the same assumption would be permissible in 
‘static’ aeroelastic problems, such as those of Section 
III (and was made in reference 3), it was not made in 
this paper because it was not necessary, a more accur- 
ate solution being possible by matrix methods without 
excessive labor. In oscillatory cases the modes may 
logically be assumed to be the natural ones of zero air- 
speed elastic vibration. In the static case, one’s best 
estimate is to use the modes of deflection under rigid- 
wing air loads, but the calculation of the rigid-wing 
loads and the corresponding modes involves nearly as 
much labor as to calculate for the elastic case directly, 
as in Section IT. 

In the present example, an airplane will be considered 
whose main wing is geometrically similar to Wing No. 1 
of the preceding Section. It will be assumed that this 
wing is extremely stiff in torsion, so that both inertia 
and aerodynamic loads due to torsional deflection are 
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negligible; this eliminates the coupling of this motion 
with the others and reduces the number of significant 
equations. It is further assumed that the wing bending 
deflection 2(s, t) is 


z= 2(s,t) = veo(s/s1) (4.1 


where v is a real constant, the amplitude of wing bend- 
ing, A is a complex constant, and ¢(s/s:) is the mode of 
fundamental bending vibration at zero air speed, de- 
fined by 


—w'*moy + (Elyss) ss = 0 (4.2 


where w is the fundamental bending frequency. 

After some labor, it will be found that the three equa- 
tions of motion—for wing bending, fuselage vertical 
motion, and fuselage pitching, respectively—take on 
the forms of Eqs. (4.3)—(4.5) below. It will be noted 
that the equation of wing bending has been multiplied 
by z and integrated, in accordance with the Rayleigh- 
Ritz theory, so that Eq. (4.3) now states the equality 
of certain integrated energy terms, which would be 
exact if ¢(s/s:) were exact for the aeroelastic case. The 
solution is being sought in the form 6 = 6e”, w = we — 
finally, the underlines are omitted. 
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Certain new notation has been introduced here, as 
defined in the Summary of Notation. The most im- 
portant is the symbol & for the operator, which, oper- 
ating on an incidence function, gives the resulting lift 
distribution, cc;. Thus @ is a kind of inverse to @; in 
fact 

a= Q@[c,] 
means the same thing as 

cc; = P[a] 
Obviously, ® depends on the choice of lift-distribution 
theory, as does @, and upon the profile, Reynolds Num- 
ber, and Mach Number. 

In Eqs. (4.3)—(4.5), ®[1] denotes the lift distribution 
for uniform incidence, @[y] and @[xs;,/b] denote the 
distributions for incidence distributions ¢(s/s,;) and 
x1,/b, and @[y’] means the distribution for an incidence 
equal to y’(s/s1). 


Eqs. (4.3)-(4.5), as written, actually apply to air- 
planes with variable m and c, although these are con- 
stant in the present example. It will be noted, how- 
ever, that the sweepback has been assumed linear, so 
that £is constant. 

It is desirable to make Eqs. (4.3)—(4.5) dimension- 
less; this is easily done by introducing the familiar 
mass ratio » and characteristic time 7 of airplane dynam- 
ics, 

(4.6) 


pw = IM/pSb, 7 = M/pSU 


and defining two dimensionless variables corresponding 
tovX andw: 


’ =dr, w' =wr (4.7) 


Discussion of Equations 


The. compatibility condition for the three resulting 
dimensionless equations is the vanishing of a deter- 
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minant involving a number of dimensionless constants 
of the wing and airplane and the dimensionless 
exponent \’. The resulting equation for \’ is a 
quartic: 


Ni+ AN3+4+ B24 OC +D=0 = (4:8) 


In working out the formulas for the coefficients A, B, 
C, and D, the following numbers have been chosen for 
the present example, in addition to the specification of 
the planform as mentioned above: 


Radius of gyration in pitch)? = //M = 0.048 } 


Z,=0 


9) 
Z./pUS = —2.52, - 


Quantities involving ® have been calculated using 
strip theory and a lift-curve slope @ = 5.05, which 
corresponds to flight at a Mach Number of about 
0.7. 

From the resulting formulas for the four coefficients, 
the following important conclusion has been drawn: 
For probable practical values of the various param- 
eters, the quartic is subject to an approximate factoriza- 


tion: 


These modifications may be important in some cases. 


42 _C¢ , 
[> +(4 5) +8] 


Pe 
[ars + —\’ + — 


B 4 = 0 (4.10) 


Moreover, the second factor in Eq. (4.10) is related to 
the customary short-period longitudinal oscillation of 
the rigid airplane, for, if w’ becomes infinite, so do B, C, 
and D, and the quartic is reduced to a quadratic: 


d’2 + (C/B)A’ + (D/B) = 0 (4.11) 


In this limiting case, as will be seen below, the ratios 
C/B and D/B depend only on the airplane stability 
parameters, and Eq. (4.11) describes the familiar 
short-period motion. 

The first factor in Eq. (4.10) represents a rapid, 
highly damped oscillation. Its frequency is slightly 
greater than the natural frequency of wing bending. 
This is, apparently, essentially wing bending vibration, 
rapidly damped by aerodynamic effects, and modified 
little by the motion of the airplane. 

The approximate factorization of the quartic means 
that there is actually no appreciable coupling of wing 
vibration with airplane oscillation but rather that the 
principal effects of elasticity are modifications of the 
short-period oscillation. 


The formulas for the coefficients appearing in the modi- 


fied short-period oscillation, except for a positive multiplicative factor, are the following: 


B = {0.0605 — 0.100» + 0.032 v? + (0.013/u) — (0.0276 v/u) + (0.0111 vw’?/u) — (M,,/pUSb) X 
[0.228 »v + (0.092 v/u) — 0.107 v2] — (15/pUSb*)[(—0.183 v/u) + (0.816/u)]} (4.12) 
C = {—0.205 + 0.0965 » — (0.049/n) + (0.0275 vw’2/n) — (My/pUSb) X 
[0.820 — 0.263 » + (0.164/u)] — (1%; /pUSb*)[1.27 — 0.976 » + (0.693/nu) + (0.228 vw’?/u)} (4.13) 
D = {0.0470 — 0.220 w + 0.110 pw — (M,,/pUSb)(0.419 + 1.26 wu — 0.098 uv + 0.228 ww’?) — 
(115/pUSb?)[—0.35 + (.576 vw’?/u)]} (4.14) 


he roots involved are, of course, 


dM’ = —(C/2B) + V(C/2B)?— D (4.15) 

As might be expected, D depends primarily on the 
static-stability derivative 17,, while C depends mostly 
on the damping in pitch, 1/;. When w’ is small—i.e., 
the wing natural frequency is low—D may become 
negative; then a positive real root appears. This 
means that the static stability has been nullified by 
wing deflection. This appears to be the same effect that 
was illustrated, in a dive pull-out, at the end of the last 
Section. The approximations of the present Section, 
including the assumption of the bending mode, are not 
accurate for this type of motion. 


On the other hand, even if static stability is main- 
tained by means of a large negative 1/,, unstable os- 
cillations may appear if C becomes negative. This 
might be described as a disappearance of damping in 
pitch due to aeroelasticity. 


Results 

In Fig. 5 are plotted typical curves illustrating the 
appearance of these undamped oscillations. For val- 
ues of » greater than about 30, it is found that w’?/y 
for undamped oscillations is practically independent 
of uw. Also, as mentioned above, these curves are 
nearly independent of  ,. 

For example, if the wing weight is 10 per cent of the 
gross weight, a critical condition may arise at 
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Fic. 5. Wing bending frequencies for divergent short-period 
oscillations. » = 30, My/pSUb = —0.047. (These curves are 
practically independent of M7.) 





M j/pSUb? = —0.085, w’?/n = 13 


The corresponding actual wing frequency is 


/ U 


a », 
wo=w/r = V13 7 hh To 
bv u 


W/Sb 


rad. per sec. (4.16) 


in consistent units, where p/p) is the density ratio. 
This critical frequency, as an example, seems sufficiently 
low to prevent trouble except possibly for extremely 


large, fast airplanes. 


Comparison with Short-Period Oscillation 


It is of interest to compare the frequency and damp- 
ing of the modified short-period oscillation of the elas- 
tic airplane with those of the corresponding rigid air- 
plane. Perhaps, although the motion is stable, it may 
be so lightly damped as to be objectionable. 

For this purpose the following airplane is considered: 


p = 30, M,/pUSb = —0.047 
y = 0.15, My pUSb? = —0.26 


Table 3 presents the results, in the form of the roots }’, 
for various natural frequencies w’, including the case 
of the rigid wing, w’ = ~, 

It appears that the damping of the motion is not seri- 
ously affected by the lowered wing stiffness. More- 
over, the frequency is reduced, so that the motion should 
be more, rather than less, comfortable to the occupants. 
In fact, there is a range of w’ for which the oscillation 
disappears and a pair of subsidences exist; this, how- 
ever, may be a condition of insufficient ‘‘static’’ sta- 
bility. 


AERONAUTICAL SCIENCES-FEBRUARY, 


1949 


Effect of Sweep 

It is also of interest to investigate the effect of sweep. 
back on this motion. To do this, at least approxi 
mately, one may simply put & and all the x(s)’s in 
Eqs. (4.3)-(4.5) equal to zero. (Obviously x, and 
X3,, cannot both be zero, but the terms involved are 
small for straight wings; the important effect is to 
eliminate coupling between wing bending and airplane 
pitching.) 

When this is done, it is found that all of the coef 
ficients are positive, regardless of w’, within the range of 
interesting airplane parameters. Hence, the possi- 
bility of either divergence or divergent oscillations is 


eliminated. 


(V) CONCLUSION 


The problem of calculating the lift distribution of an 
elastic airfoil with sweep has been set up and its solu- 
tion by matrix methods indicated. It appears that in 
many practical cases it will be unnecessary to assume 
the modes of elastic deflection. The choice of the aero 
dynamic theory, as well as the exact nature of the elas 
tic theory employed, has been left free, provided only 
that these are linear and one-dimensional; thus, the 
elastic properties may be expressed by linear one-di 
mensional influence functions. 

Sample calculations of lift distribution, torsional 
divergence, and aileron effectiveness by this method 
have confirmed the known behaviors of sweptback and 
sweptforward wings. In addition, a case of dive pull- 
out at constant speed has illustrated a more compli- 
cated case where several variables are involved, which 
are related by the airplane's dynamical equations. It 
appears that the destabilizing effect of elastic deflection 
under air load is partly counteracted by the effect of 
deflection under inertia load. This argues for increased 
weight distributions for wings. 

On the other hand, an investigation of a simplified 
dynamic case involving wing bending and airplane 
longitudinal motion, with sweepback, shows possibili- 
ties not only of divergence but also unstable oscillation 
when the natural frequency in bending is low. In this 
calculation it was practically necessary to assume the 
wing bending mode to simplify the work. Both insta- 
bilities are related to the sweepback, for they do not 
seem to occur with a straight wing. 

Throughout the paper it has been necessary to resort 
to numerical examples to illustrate the methods pro- 
posed. Naturally, the conclusions drawn from such 
examples may not be general and should be regarded 


TABLE 3 
Characteristics of Modified Short-Period Oscillation 
, ‘ 





@ A 

o0 —3.92 = 5.221 
100 —3.74 + 4.397 
55 —3.35 = 2.302 
42 —0.586, —5.65 
39.5 0, —5.9 
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accordingly. The methods themselves are rather gen- 
eral, however, and will handle problems considerably 
more complex than the rather simplified numerical 
cases treated here. 


Appendix 


He AERODYNAMIC MATRIX [@] IN VARIOUS APPROXI- 
MATIONS 
The matrix [@] is defined by the equation 
ta} = [@] {cr} (A.1) 
A) Strip Theory 

a(s) = c;,(s)/a(s) (A.2) 
Hence, [@] is the diagonal matrix whose members are 

the values 1/a(s,) (pb = 0, 1,2,...7). 


B) Prandtl Theory 
The Prandtl integrodifferential equation, in trigono- 
metric-series form, is 


n 


Lb ’ 
¢, = A, sin mdb = 
‘ 


m=1 
n 


j ) sin mol 
asa — mA» — 
] sin 3 \ 


(A.3) 
m=1 / 
where ? = cos~!(2y/d). 

Glauert’s solution"! consists of writing this equation 
out for 2 values of y, using known numerical values of 
¢,a, and a for these values of y and then solving the 
set of linear equations for A;,...A,. In matrix nota- 
tion, this process can be expressed as follows: By 
Fourier’s theorem, 


b/2 
l , , 
A, = -—— cc, sin md(dd/dy)dy 
2rbJ/ -0/2 


(a) « I 
Ai 3G 


or 


[=] [S] [ce] [1/sin #] {c;} (A.4) 
where [=] denotes the square matrix whose members 
are the values of Z,, = sin(md,); m, p = 0,1, 2,...m; 
here 3, denotes cos~'(2y,/b), where y, is the value of 
Also [1/sin #] 
tepresents the diagonal matrix whose members are 
l/sin 3. 


Now, according to Eq. (A.3), 


[a]—"{c,} + [1/sin 8] [E]"[n] {A} 


f 
ja} = 


ll 


((a] + 5 f1/sin OE} in LEIESIe] X 


[1/sin #]){c,} (A.5) 
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where [a]—! is the reciprocal of [a], thus equal to [ @] of 
the Strip Theory, above, and [=]? denotes the trans- 
posed matrix of [=]—i.e., the same matrix but with 
rows and columns interchanged. The diagonal matrix 
[7] has the values 1, 2, 3, . . . m along its diagonal. 

Hence, the matrix [@] for this case is the entire ma- 
trix enclosed in parentheses in Eq. (A.5). 


(C) Weissinger Approximation 

Weissinger’s approximation consists in adjusting the 
circulation distribution to make the slopes of the stream- 
lines at the three-quarter-chord points compatible with 
the given angle-of-attack distribution. The downwash 
at this location is calculated by integration over the 
sweptback bound vortex and its wake. The result can 
be written 


—. 1 - "b/2 ~ \ J 1 1 ) 
a= ee (n S - yt gh n) dn (A.6) 


where I is the circulation, equal to '/2U/cc,;, and L(y, n) 
is a function calculated by Weissinger,* which depends 
on the wing geometry. 
If, again, is expressed as 
n 


2Ub ¥* A,, sin md 


aw 
1 


the first term in the integral in Eq. (A.6) takes on the 
same form as the second term in Eq. (A.5). The 
second term in Eq. (A.6) can be expressed as follows: 


b/2 *b /2 
7 l 7 
“ ff I'’(n) L(y, n)dn = — / cc, Ldn (A.7) 
4nrU J -0/2 : SrJ —b/2 


where Ly denotes OL(y, n)/On. Thus, 


[@] = (( ‘sin 8) [=]? [n][E][1/sin 8] — 


l 
mb 
b 
Us} elt] (A.8) 


Here [L,] is a square matrix whose terms are the values 
Li(¥p, 1); P»r =0,1,...0. 
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* See reference 7, Eqs. (36) and (37). L(y, 7) is equal to 
L¢(¥, 7) multiplied by the aspect ratio \ and divided by the half 
span, }/2. 

(Continued on page 119) 








A Simple Approximate Formula for 
Effective End-Fixity of Columns 


N. M. NEWMARK* 
University of Illinois 


SUMMARY 


This paper gives simple approximate formulas for numerical 
values of the end-fixity coefficient of a uniform bar, with elastic 
restraints at the ends, subject to compressive forces. These 
formulas were developed as part of a project sponsored at the 
University of Illinois by the Consolidated Vultee Aircraft Cor- 
poration and were first reported on March 15, 1944. Apprecia- 
tion is expressed to the Consolidated Vultee Aircraft Corporation 
for permission to publish the formulas in this paper. 


DISCUSSION 


| PR yerenpae A BAR of span L, modulus of elasticity E, 
and uniform moment of inertia J, subjected to an 
end thrust of magnitude P. The bar is considered to 
be initially straight. The ends of the bar are subject 
to elastic restraints against rotation either by springs or 
by continuity with other bars; however, the ends of 
the bar are not free to deflect. 

In the following, the subscript 1 refers to the left end 
of the bar and the subscript 2 to the right end of the bar. 
The stiffness against rotation of the restraint at the left 
end of the bar is characterized by the equation: 


M, = Ridy (1) 


where J, is the moment exerted on the bar by the 
restraint and ¢; is the rotation at the left end of the 
bar. At the right end of the bar a similar relationship 
holds: 


Mz = kode (2) 

The critical end load P can be stated as a coefficient 

of end-fixity, C, multiplied by the ordinary Euler load 
for pinned ends: 

P = C(r*EI/L?) (3) 
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It is convenient to use dimensionless quantities m and 
n related to the stiffness k; and k in the following way 


m = k,(L/E), n. = k,(L/EI) (4 


Exact formulas relating m, m2, and C can be written, 
but these are relatively complicated and do not give a 
simple picture of the mechanism of the behavior of the 
bar. Asa result of numerical calculations from the ex- 
act formulas, and with simplification of these formulas, 
approximate equations have been derived and are 
given here. 


The most general formula is given in Eq. (5a): 


Pos (SS \(E+ ss) (5a) 

mw? + 2,/\mr? + 2n2 
This equation is in error by less than 4 per cent and is 
valid for infinite as well as for finite values of m; and m, 
The values of C given by this equation are generally 
less than the true values with the maximum errors oc- 


curring in the range of (#7; + m2) from 10 to 20. For 
symmetrical restraints, 2; = m2, the error is always less 





than 2!/, per cent. 
Almost as good an approximation is given by the fol- 
lowing equation: 


Ce (ese Co ov) (5b) 
1 + 0.2m,/\1 + 0.2n2 


which may be simpler to use for small values of n. 
When (m + mz) is less than 1.0, Eq. (5a) or (ob) 
may be simplified with the results shown in Eq. (6): 
2 
C2 1+ —(m + m) =14+0.2(m + m) (6) 
gt 
Eq. (6) is in error by less than 2'/, per cent when (m + 
nz) is less than 1.0. 
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Notes on Column End-Fixity Coefficients 


G. C. BEST* 


Consolidated Vultee 


ABSTRACT 


The theory relating elastic end restraint to end-fixity coefficients 
for beams of constant section is developed in a somewhat dif- 
ferent manner than heretofore. A new and simple graphical re- 
lation results from which arises an approximate formula of an 
accuracy adequate for engineering purposes. Also an approxi- 
mate formula discovered by N. M. Newmark is shown. 


INTRODUCTION 


| Bipagner-eone LITERATURE has been presented on 
the subject of column end-fixity coefficients.'~’ 
The formula relating such coefficients to the degree of 
elastic end restraint has been derived in various pre- 
vious papers.*» * In the following analysis the slope 
deflection method is used to obtain this relation. As 
is customary, it is assumed that the Euler bending 
equation applies and that deflections are sufficiently 
small to permit the substitution of dy?/dx? for M/EI. 
Also, all strains are assumed within the proportional 
limit. The column considered is of constant section. 


Basic THEORY 


Assume a column loaded as shown in Fig. 1. The 
elastic restraint at end 1 is such that a rotation of ; 
radians causes a moment of J/,’ in.Ibs. Similarly, 
at end 2, or: 


M,’ =s Kh, M,’ a Kobe (1) 


Defining a positive moment as one that tends to turn 
a joint clockwise and a positive rotation as a clockwise 
rotation, the slope-deflection equations for the column 
of Fig. 1 become 


0 = (Ki+ S’)@, + pace 


0 = S'Ch + (Kz + 5S’), (2) 


where 
S’ = beam stiffness of column with other end fixed, 
including effect of axial load 
C = carry-over factor, including effect of axial load 


For Eqs. (2) to have solutions other than zero, the 
determinant must be zero; hence, 


(Ky + S’) (Ke + 5’) = S”C? (3) 


It is possible to write S’ = oS, where S = beam stiff- 
ness, other end fixed, when no axial load is applied 
(iie., S = 4EI/L), and the factor o constitutes a cor- 
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rection factor taking the effect of the axial load into 


account. In the notation of reference 1: 
¢ = 38/(48? — a*)| : 
C = «/28 f (3a) 
where 
a = 6[(L/j) cosec (L/f) — ILD gy, 


3 = 3[1 — (L/f) cot (L/)V/(L/i)* J 
and m/c = L/j where c = end-fixity coefficient. On 
this basis Eq. (3) becomes, dividing through by S° 

[(Ki/S) + o][(K2/S) + ¢) = °C (4 
It is now advantageous to introduce the terms 


dy = § ‘(Ky + 5), dg = S/(Ke a S) (5) 


where d; and d, are the distribution factors to the col- 
umn at ends 1 and 2, respectively, which would exist 
if the column were not under axial load. From Eq. 
(5) it is possible to write 


Ki/S=(1—4d)/d, Ks/S = (1 — ds)/d2 (6) 


Substituting in Eq. (4) and letting A = 1 — o and 
h? = o°C? — A’, for brevity, gives: 
1 — A(d; + dz) = deh? (7) 
To facilitate plotting let 
d = (d,; + d)/2, ¢ = d,d2/d* (8) 
Eq. (7) then becomes 
1 — 2Ad = ed*h? (9) 


Solving Eq. (9) for d and neglecting the false root 
gives 


—A+ Va? + ch? 














d = = 
eh? 
—A+ V(1— A? + °C (10) 
eh? 
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The last indicated form shows that the expression under 
the radical will never become negative, since 0 < € < 1. 


SPECIAL CASES 


When ¢ = | the elastic restraint at both ends is the 
same. For this case Eq. (4) becomes 
(K/S) + a= 0C 
or 
[(1 — d)/d|+ea=0C 
or 


d = 1/(A + oC) (11) 


When ¢ = 0, either d; or d; must equal zero—i.e., one 


end must be fixed. For this case Eq. (9) yields 
a= T/24 (12) 


0 with- 
out assuming negative elastic restraint at the end that 


For d > 0.5 it becomes impossible to have e = 


For example, supposing end 2 fixed, then 


d=d, = 


is not fixed. 
(1/2) [S/(A, + S)] 
which cannot exceed 0.5 unless A, < 0. 


So as not to interrupt continuity, both Eqs. (11) 
and (12) are plotted in Fig. 2, though negative end re- 
It is 
possible to interpolate linearly for «€ between these 


straints are not considered to actually occur. 
curves with an accuracy (errors <0.15 per cent) greater 
than that to which the graph can be read. 


APPROXIMATE FORMULAS 


The appearance of the graph of Eq. (11) suggest the 
use of a parabolic approximation neglecting «. Such 


a formula, giving greatest errors of <2 per cent, is: 
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Vc = 2 — [d(6 — d)/5] (13 


Another approximate formula, for c rather than Vc, de 
veloped by N. M. Newmark, is 


¢ = (Ai + 0.4)(A2 + 0.4)/(A1 + 0.2)(Ag + 0.2) (14 


where A; = S/4K, and dy = S/4Ko, 

The greatest errors occurring from the use of this for 
mula are <4 per cent which, since it is for c, instead of 
Vc, makes it of comparable accuracy to Eq. (13). The 
Newmark formula has the advantage of symmetry 
Also, if the end stiffness at end 1, say, is measured by 
the end-fixity coefficient (c,;) that would result if both 
ends were of like elastic restraint and similarly for end 
2, Eq. (14) yields the simple result 


/ 


c= ¥ C\C2 (15 


ILLUSTRATIVE EXAMPLE 


The application of the formulas is illustrated by the 
following problem. In Fig. 3 is shown a column with 
end restraint. It is assumed to buckle in the plane of 
the paper. The relative of the 
members, with the axial load on (1—2) ignored, are: 


stiffnesses various 


Member Stiffness 


1-3 49/50 = 0.8 
1-2 10/30 = 0.33% 
9-4 100/50 = 2.0 


The distribution factor at end 1 on this basis is 


d, = 0.333/(0.8 + 0.333) = 0.292 


and at end 2 is 


Hence 


and 
e = (0.292)(0.143)/(0.218)? = 0.88 


. - a / — 
From Fig. 2 one obtains Yc = 1.755. Hence, 


P = cx*(El/L*) = (1.755)*[a*10*/(30)? = 338 Ibs. 
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COLUMN END-FIXITY 


The first approximate formula gives 


Vc = 2 — (0.218/5)(6 — 0.218) = 1.748 


and the Newmark formula gives 


\1 = 0.333/(4)(0.8) = 0.104 
hy = 0.333/(4)(2.0) = 0.407 
or 
c = (0.504) (0.4407) /(0.304) (0.2407) = 3.04 
or _ 
a/c = 1.75 
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Analvsis of Eccentrically Loaded Riveted 
or Bolted Joints 


G. M. MONROE* 
United States Navy 


SUMMARY 

The semigraphical analysis developed here is applicable to all 
riveted or bolted joints not stressed beyond the elastic limit. 
The stress on each rivet or bolt is obtained by plotting one point 
on a scale drawing of the joint and multiplying the distance from 
this point to each rivet or bolt by a constant factor. The co- 
ordinates of the point and the multiplying factor depend only on 
the loads on the joint, the total rivet area, andthe polar moment 
of inertia of the rivet areas. The application of this method is 


simple and brief. 


INTRODUCTION 


A SHORT AND SIMPLE METHOD of determining the 
rivet or bolt loads in eccentrically loaded riveted 
or bolted joints is presented herein. This method has 
distinct advantages over the usual method as described 
in most handbooks and textbooks. It requires a mini- 
muin of computation, which leads to a graphical solu- 
tion that can be performed on any sheet of coordinate 
paper, or even on the original working drawing, without 
the use of drawing instruments or protractors. 

The solution as usually presented requires determina- 
tion of the radial distance and direction from the center 
of resistance of the joint to each rivet of the joint. Us- 
ing this distance, the loads on each of the rivets must 
be computed analytically, the results occurring as two 
normal components along the two chosen axes and a 
third component acting normal to the radial direction. 
The total rivet load must then be found by combining 
these three force vectors, a tedious operation that must 
be performed individually for each rivet in the joint. 

The graphical method will be developed and ex- 
plained in general terms, using Fig. 1 to indicate nomen- 
clature. The joint of Fig. 1 will be analyzed to illus- 
trate the method advocated. & 

The center of resistance of a riveted joint is the cen- 
troid of the rivet areas. In Fig. 2 the applied loads 
have been replaced by a new force system acting at the 


center of resistance. 


SYMBOLS 


P,;P, = applied forces at center of resistance 

M = applied moment at center of resistance 

F = load on rivet, lbs. 

A = area of cross section of rivet (4; = total area of 
rivets) 

Ss = shear stress on rivet 
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R = radial distance of rivet from center of resistance 
X, Y = coordinates of rivet measured from center of re. 
sistance 
a,b,etc. = subscript, denotes rivet 
Prime = stress due to applied forces 
Double 
prime = stress due to applied moment 
Ip = polar moment of inertia of rivet areas about center 
of resistance 
Iz, ly = moment of inertia of rivet areas about X or Y 
DISCUSSION 


Basic formulas are derived as follows: 

The components of the load, P, and P, acting through 
the center of resistance, cause equal stress on all rivets. 
The stress components for any rivet are given by 

== 2 . Pee / 
S,’ = P,/A,; S,’ = P,/A;s (1 
The stresses caused by J/ acting about the center of 
resistance are obtained as follows, making the usual 
assumption that the stress on each rivet varies as its 


distance from the center of resistance. From this pro- 
portionality it is seen that 


Sa"/Ra = Sy”"/R, = S,."/R, = ete. (2) 


The applied moment is equal to the sum of the moments 
of the various rivets—that is, 


M = S,"AgRq + Sp"AsRo + ... + Sn"AmRn (3) 
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From Eq. (2): 
Sa" Ro, 


YW 
Sp =e 


, 
Ra 


Sa’Re. 
 ” 


S,” = 5S,” = etc. 


and if these values are inserted in Eq. (3), there results: 


S,"i 5,"Re 

R, AR +...t+ — 
(Sq”/Rq)(AgRq? + ApRy? + ... + AmRm?) 
(Sq”/Ra) (Io) 





M = eige + AnRmn 


Solving for the stress on a rivet due to the moment M, 
S” = MR/Iy (4) 


The direction of S” is normal to the radial distance, and 
from Fig. 3, by similar triangles, the components of S” 
are: 


S," = (Y/R)S"” = (M/Ip) Y 
Sy” = —(X/R)S" = —(M/I)X (5) 


By adding the stress from Eqs. (1) and (5), the total X 
and } stress components are found. These components 
are given by the following basic equations. 


Se = (P,/A,) + (M/Io) Y (6) 
S, = (P;/Ad) — (M/Lo)X (7) 
The terms in the parentheses are constant for a particu- 


lar problem and may easily be determined. Note then 
that S, is a function of only Y, and S, is a function of 


i 


only Y. For graphical application Eqs. (6) and (7) are 
rewritten as follows: 
S-/(M/Io) = [(P2/M)U0'Av] + ¥ (8) 
Sy/(M/Io) = [(Py/M)(0/A1)] — X (9) 


Eqs. (8) and (9) are general for any riveted or bolted 
joint and, as given above, are for a joint whose rivets 
are in single shear. If the rivets are in double shear, 
the stress will be one-half of that indicated by the equa- 
tions. 

The complete analysis is performed as follows: If 
the Y and Y coordinates of all the rivets are plotted on 
a sheet of coordinate paper, using the center of resist- 
ance as the origin, it is possible to locate a point Q such 
that the distance from Q to any rivet is proportional to 
the stress acting on that rivet. This follows from the 
linear relations shown in Eqs. (8) and (9). The point 
Q must be so located that its XY coordinate is P,Jo/MA,. 
Then the XY distance (parallel to the Y axis) from Q to 
any rivet is (P,J)/MA,) — X, which from Eq. (9) is 
S,/(M1/Io) for that rivet. Similarly, the ])’ coordinate 
of Q must be —(P,J)/MA,). Then the J] distance 
(parallel to the Y axis) from Q to any rivet is —(P;lo + 
MA,) — Y, which from Eq. (8) is —S,/(\//Jo) for that 
rivet. (See Fig. 4.) The distance measured directly 
from Q to a rivet is proportional to the total stress on the 
rivet. That is, 


-_ ( S, \'+ (Ge) 
w/t \ar/r W/To 


where SS is the total stress on the rivet. The necessary 
steps in locating Q and determining the stress are out- 
lined below. The construction is indicated in Fig. 4. 
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(1) Plot the X and Y coordinates of each rivet on co- 
ordinate paper, using the center of resistance as the 
origin. 

(2) Compute the coordinates of Q as indicated above 
and plot Q to the same scale as that used in Eq. (1). 

(3) Measure the distance from Q to each rivet and 
multiply by the factor .1//Jo if the rivets are in single 
shear, J//2/) if in double shear, to obtain the total 
stress on the rivet. The direction of the stress is normal 
to the line connecting the point Q and the rivet. 

(4) If the force on the rivet is desired, the stress 
obtained in step (3) should be multiplied by the area 
in shear of the rivet. 

These four simple steps are sufficient to solve any 
riveted or bolted joint. If a scale working drawing ot 
the fitting is available, it is not necessary to replot the 
rivets on coordinate paper. The point Q may instead 
be located on the drawing, using the scale of the drawing 
to plot and measure distances. The most highly 
stressed rivet is the rivet that is the greatest distance 
from Q. Therefore, if only the margin of safety of the 
rivets in shear is desired, only one distance need be 
measured, 

The following analysis of the joint shown in Fig. | is 
included to illustrate the use of this method. 

Total Rivet Area: 

A of '/sin. rivet is 0.0491 sq.in.; A of */16-in. 
rivet is 0.0276 sq.in. 
A, (total rivet area) 8(0.0491) + 4(0.0276) = 
0.5032 sq.in. 
Center of Resistance, <, ¥: 
& = (1/0.5032) [38(0.0982) + 2(0.0982) + 
1(0.0982)] = 1.172 in. 
y = 0, by symmetry. 
A pplied Loads: 

P. = 19,800 Ibs.; 
Applied Loads Transferred to Center of Resistance: 
1,000(1.625 + 1.172) — 


P, = 


u 


~4,000 Ibs. 


AJ = moment change = 
19,800(0.125). 

M = 9,120 in.Ibs. (clockwise). 

Polar Moment of Inertia: 

I, = 8(0.0491) (2.5)? — 2(0.0276) (1.5)? —2(0.0276) X 
(0.5)? = 2.59 in.‘ 

I, = 2(0.0491)(32 + 2? + 
0.680 in.* 

Io = I, + I, = 3.27 in. 


1°) — 0.5028(1.172)? = 
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Computation of Constants: 
P,Ip/MA, = —14.12 in.; P,lo/MA, = —2.86; 


M/Iog = 2,790 lbs. per sq.in. per in. 


The point Q in Fig. 5 is located using these constants 
The distance from Q to each rivet is measured and re- 
corded in the third column of Table 1. The shear stress 
and load on each rivet are computed in the fourth and 


fifth columns. 


TABLE 1 


Rivet S- S, Lbs. SA = 

Rivet Area M/TIo per Sq.In. F 
, 0.0491 16.55 46,200 2,270 
b 0.0491 16.75 46,800 2,300 
c 0.0491 16.90 7,200 2,320 
d 0.0491 17.10 47,700 2,340 
e 0.0276 16.30 45,500 1,256 
f 0.0276 15.16 42,300 1,168 
g 0.0276 14.18 39,600 1,093 
h 0.0276 13.26 37,000 1,021 
i 0.0491 [3.23 34,200 1,680 
j 0.0491 12.00 33,500 1,645 
k 0.0491 11.78 32,800 1,610 
m 0.0491 11.68 32,600 1,600 
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of the Libraries of the 


Institute of the Aeronautical Sciences 





The services of the Libraries are available to all mem- 
bers of the Institute, to Corporate Members, to advertisers 
in the AgrRoNnauTICAL ENGINEERING Revisw and AERO- 
NAUTICAL ENGINEERING CaTALoa, and, under usual library 
limitations, to the public. Four specialized services are 


available. 


The Paul Kollsman Lending Library 


This lending library service makes available, without 


| charge, the latest and more important aeronautical books. 


Members may request the loan of any aeronautical or 
technical book they wish to borrow. Through an exchange 
agreement with the Engineering Societies Library, any 


| book on general engineering may be borrowed from its great 
| collection of over 160,000 volumes. 


A photostating service is available at usual library rates. 
Applications for membership in the library and further 
information will be sent on request. 


The W. A. M. Burden Reference Library 


This reference library contains over 12,000 aeronautical 
books, magazines, pamphlets, and reports gathered from 
world-wide sources and is one of the most complete aero- 
nautical libraries in the world. Material from this library 


| is not available for loan but may be used for reference pur- 


| poses. 


The Pacific Aeronautical Library 


6715 Hollywood Boulevard 
Los Angeles 28, California 


Established in cooperation with the aircraft companies 
the library serves. The leading aircraft companies in or 
near Los Angeles participate in its support and operation. 

This service library for aeronautical research is available 
to the public for reading privileges. Source material in- 
cludes aerodynamic and structural research reports, as well 


as books on drafting, production methods, history, and al- 
lied sciences. It furnishes books, periodicals, and pamphlet 
material to the participating aircraft companies to supple- 
ment their engineering libraries. 


Technical Information Service 


This service has experienced personnel under the super- 
vision of trained aeronautical engineers to compile any in- 
formation desired. The services range from listing special- 
ized reference books to the preparation of exhaustive 
bibliographies, digesting of reports, and general surveys of 
any aeronautical subject. Some of the available services 
are: 


Bibliographies on any aeronautical subject. 

Reports on any aeronautical subject. 

Digests of aeronautical books, papers, periodicals, and refer- 
ences. 

Translations. 

Engineering investigations of special aeronautical subjects. 

Biographies of individuals engaged in aeronautics. 

Photostats of any aeronautical or general engineering 
material. 

Microfilms made on special order. 

Photographs made from the Institute’s photographic collec- 
tion. 

Drawings and tracings made. 


In addition to the services mentioned any commission 
which comes within the scope of the Service will be ac- 
cepted. Special arrangements may be made for work re- 
quiring several weeks or months. 

Translators are available for accurate transcriptions of all 
foreign language data. Translations are carefully edited 
by trained engineers. 

Reproductions of any material in the collections of 
the Institute may be ordered at standard photostat rates. 
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Stringer Panel and Multicell Wing 
Construction 


Kenneth P. Buchert 
Langley Aeronautical Laboratory, N.A.C.A., Langley Field, Va. 
September 24, 1948 


I‘ THE PAPER ‘Efficient Application of Stringer Panel and 
Multicell Wing Construction’? (JouRNAL OF THE AERONAU- 
TICAL ScrENcES, Vol. 16, January, 1949), Mr. Gerard has de- 
scribed a method for approximately locating the boundary be- 
tween the region in which stiffened panel construction will more 
efficiently meet the design conditions and the region in which 
multicell construction will be more efficient. Because the anal- 
ysis necessarily does contain important assumptions and ap- 
proximations, perhaps some additions are in order to Gerard's 
already well-considered discussion of the possible shifts in the 
boundary due to these approximations. 

First let us consider those assumptions and approximations 
which favor the stiffened panel. 

(1) Restriction of Analysis to One Type of Multicell Construc- 
tion.—There is available a vast backlog of testing and experience 
with many different types of stiffened panel construction which 
makes possible the choice of an optimum stiffener configuration. 
This backlog has as yet no counterpart in multicell construction. 
Already designers are proposing different types of multicell con- 
struction employing for example, corrugated shear webs, shear 
webs that act as diagonal tension beams, single stringers between 
webs, or sandwich skin between webs. After these types have 
been evaluated, a comparison can be made between stiffened 
panel construction and the optimum type of multicell construc- 
tion, rather than the particular type assumed by Mr. Gerard. 

(2) Negligibility of Rib Weight in Stiffened-Panel Design.— 
Authorities differ regarding the desirable relative weight of ribs 
and stiffened panel. kK. J. Farrar, in a Bristol Aeroplane Com- 
pany paper, came to the conclusion that, for optimum design, 
rib weight should be approximately one-half the panel weight. 
Inasmuch as the relative amount of rib weight required may be 
expected to increase as the depth of the wing increases, Farrar’s 
rule of thumb probably should be considered most nearly ap- 
plicable at high values of the ratio b,,./nb, in Figs. 5 and 6 of Ger- 
ard’s paper; or, in other words, if rib weight were taken into con- 
sideration, the stringer panel would appear somewhat less favor- 
ible for thick wings than indicated. 

(3) Negligibility of Stiffener Height—As Gerard points out, if 
the height of the stiffeners is completely overlooked, a point is 
reached at which the stiffeners on the upper and lower skin will 


touch. The reduction in effectiveness of the stiffened panel to 


resist bending of the wing, due to the fact that the center of 
gravity of the panel cross section lies nearer the neutral axis of 
the wing than the center of gravity of the skin alone, is such that 
the stiffened panel would appear somewhat less favorable for thin 
wings with wide rib spacings than indicated by Fig. 6. 
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There are, on the other hand, assumptions in the analysis which 
favor the multicell type of construction. These are: 

(1) Negligibility of Shear and Crushing on the Web in Multicel] 
Construction.—If the webs are considered to be subjected to ver- 
tical shear and to crushing because of the spanwise bending of 
the wing, the bending buckling stress of the web is reduced. If 
this reduction were taken into account the multicell construc. 
tion would appear somewhat less favorable. 

(2) Negligtbility of Side Support on the Strength of the Stiffened 
Panel.—This is the common assumption to make, and in general 
it is a fair one because the panel width is usually so great com- 
pared to the panel length, and the transverse stiffness so small, 
that the effect of side support is negligible. For cases such as the 
two in Fig. 6, however, which are the most favorable to the multi- 
cell, effect of side support may not be negligible, since the panel is 
approximately square (mb; = 25, L = 24o0r30). In other words, 
toward the lower right-hand corner of Fig. 6, evaluation of the 
effect of side support might make the multicell construction ap- 
pear somewhat less favorable. 

(3) That the Multicell Is an Integral Monolithic Structure. — 
The N.A.C.A. investigation of the effect of riveting on the 
strength of stiffened panels indicates that, in order to approach 
the strength of an integral structure, extremely strong rivets 
closely spaced are required. As yet it has not been demonstrated 
that any combination of rivet diameter and pitch produces the 
same local buckling strength that is achieved by an integral 
structure. This fact, together with the lack of experimental 
verification of the strength of multicell construction, suggests 
that perhaps the multiceil is somewhat less favorable through- 
out the range of Figs. 5 and 6. 

The foregoing comments are not intended to be criticisms of 
Mr. Gerard’s paper. It would have been impossible to include 
all the factors mentioned in the pioneering comparison he has 
made, 





Rotary Compressors 


Hunt Davis,* D. G. Traver,f and A. M. G. Moodyt 
Elliott Company, Jeannette, Pa. 
November 19, 1948 


i THE PAPER ‘‘The Relative Merits of Rotary Compressors” 

(JOURNAL OF THE AERONAUTICAL SCIENCES, Vol. 15, October, 
1948), the authors have given serious thought and study to some 
fundamental problems of compressors. They have, however, 
overlooked several important points that tend to vitiate their con- 
clusions. 

In the first place, the concept of an efficiency envelope must be 
looked at much more closely than the authors have done. They 
use the term rotary to include two essentially different types of 
compressors—i.e., the dynamic steady flow type and the positive 
displacement nonsteady flow type. Common usage of the term 
“rotary” as applied to compressors reserves it for rotating posi- 
tive displacement compressors. Unfortunately, there is no com- 
monly used simple term to cover the dynamic type. This may 
account in part for the authors’ failure to distinguish between 
them. 

It will be readily realized that the mechanism of compression 
in the two types is fundamentally different. At one point the 
authors say: ‘It is reasonable to suppose, therefore, that a ma- 
chine of comparatively high loading will be inherently less effi- 
cient than a machine producing a smaller head per stage, unless 


* Division Engineer, Aerodynamics Division. 
t Research Engineer, Analytical Section. 
t Division Engineer, Blower Division. 
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READERS’ 


the path of flow through the latter is considerably more com- 
plicated.” In the case of a positive displacement compressor, 
however, the term “‘flow path” is hardly applicable. It might be 
said that the equation of continuity and the Euler and Bernoulli 
equations do not apply. There is, therefore, no reason to expect 
a continuous curve of maximum efficiency versus load factor cov- 
ering compressors both of the positive displacement type and the 
steady flow type. 

Within one general type (i.e., ““dynamic’’) an efficiency en- 
yelope undoubtedly exists. since all subtypes (centrifugal, mixed- 
flow, axial-flow, propeller) can be produced by continuously vary- 
ing the geometry. 

There is, however, no justification for the statement made in 
the conclusions to the effect that ‘‘the maximum efficiency . . . in- 
creases from one type of compressor to another of lower char- 
acteristic load factor.’’ This conclusion, of course, is related to 
the authors’ thoughts about flow path. If tLis conclusion were 
correct, the propeller type of fan would be more efficient than 
any other type. According to the authors’ own figures, this is 
not thecase. After all, efficiency is a relation in which both losses 
and useful work are involved, and the mere fact that losses may 
be low in a particular type of compressor does not mean that 
efficiency will be high. 

The data of Fig. 2 have enough spread in Reynolds Numbers 
to warrant some adjustment for the Reynolds Number effect. 
Experimental data on modern centrifugal compressors!:? indi- 
cate that the variation in small stage efficiency with load factor 
should show a peak of about 85 per cent for a load factor of 5, 
falling off to 78 per cent at a load factor of 8. These values are 
for Reynolds Number 150,000 based on hydraulic diameter of the 
impeller passage. Campbell and Talbert’s data! are typical 
fora Reynolds Number considerably lower. The envelope shown 
by the authors seems to fall off too rapidly as load factor increases. 

The data on Lysholm compressors are definitely faulty. Pub- 
lished data are available? showing efficiencies in the neighbor- 
hood of 80 per cent for load factors ranging between 20 and 40, 
Thus, it is apparent that the Lysholm compressor possesses char- 
acteristics that compare favorably with axial-flow and centrif- 
ugal compressors for many applications including the gas-tur- 
bine power plant. 

The authors confine themselves to a consideration of efficiency. 
It must not be forgotten, however, that other characteristics, 
particularly stability, may be of major importance in a particular 
compressor application and that stability may influence efficiency. 
Axial-flow compressors are so limited in their stable flow range 
that their usefulness for industrial applications is seriously im- 
paired. Centrifugal compressors are so widely used in applica- 
tions where stability is important that in a sense they are not 
“bred” for maximum efficiency. 

With regard to Eqs. (1) and (2), it would seem that they should 
include, as independent variables, y, the ratio of the specific heats, 
as well as the dimensionless ratios describing geometrical propor- 
tions. The geometrical ratios will serve to show the dependence 
of ¥ and y upon the geometrical proportions when all the other 
variables are held constant. 

In connection with the discussion of Reynolds Number and 
referring to reference 7 in the paper, it will be found that the Reyn- 
olds Number is based on the peripheral velocity of the blade and 
not the mean axial velocity of the gas. Fig. 1 in the paper is 
taken directly from reference 7 and the correct definition of 
Reynolds Number applies. 

With regard to surface roughness, it has been our experience 
that boundary layers are almost always so thick in centrifugal ma- 


Campbell, Kenneth, and Talbert, john E., Some Advantages and Lim#- 
lations of Centrifugal and Axial Flow Aircraft Compressors, SAE Journal, 
Trans., Vol. 53, No. 10, p. 607ff, October, 1945. 

*Karrer, W., The Experimental Gas Turbine at the Ocerlikon Works in 
Switzerland, The Engineers Digest, Vol. 5, No. 5, p. 229ff, July, 1948. 

* Wilson, W. A., and Crocker, J. W., Fundamentals of the Elliott-Lysholm 

ompressor, Mechanical Engineering, Vol. 68, p. 514, 1946. 
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chinery that the effect of roughness is obscured. This means 
that the friction losses in centrifugal compressors vary inversely 
with the fifth root of the Reynolds Number, approximately. 
This agrees with the loss analyses done some time ago on hydrau- 
lic pumps and turbines. 

The Reynolds Number curve reproduced from reference 7 of 
the subject paper not only shows the Reynolds effect but includes 
an aspect ratio effect as well. There appears in Fig. 27 of refer- 
ence 7 an illustration of the group of axial-flow compressors from 
which these test results were obtained. They appear to have 
identical blade lengths of about 1 in. and constant diameter ratios. 
Consequently, these compressors are not geometrically similar. 
In these tests, the chords varied from 6.50 to 0.65 in. at the outer 
diameter, resulting in a large variation in aspect ratio. 

There is an error in the denominator of Eq. (3). 
the term (o;/pi:) there should be no division sign. 

To sum up, the writers would submit the following conclusions 
to repiace the third and fourth conclusions of the authors’ paper: 

(1) Within the range of the dynamic steady flow types of com- 
pressor there undoubtedly is a continuous efficiency envelope for 
a fixed value of Re, although this curve may not have a single 
maximum point. 

(2) On the basis of recorded data the centrifugal compressor 
has approximately 2 to 3 per cent less maximum efficiency than 
the axial-flow compressor of comparable size. 


Following 





Errata 


Chi-Teh Wang 

Associate Professor of Aeronautical Engineering, New York Uni- 
versity, New York, N.Y. 

November 17, 1948 


N THE WRITER’S PAPER entitled “‘Variational Method in the 
Theory of Compressible Fluid” (JouRNAL OF AERONAUTICAL 
Sciences, Vol. 15, No. 11, pp. 675-685, November, 1948) there 
are some errors and misprints that should be corrected. They 
are as follows: 
Page 681, Right Column: 
Line6: “y = 30ry = 4” should be “y = 1.50r y = 1.333.” 
Line 11: ‘‘taking four constants A3,...” should be “‘taking 
six constants Ay, Ais, Aa, Ass, Aus, An.” 


Page 682, Right Column: 


Line 4: ‘Ay’ = —0.03727” should be ‘‘A33’ = —0.03727.” 
. “ob ” “od 1 
Line 6: - = should read: = 0.417 + cos 6 + 
do do r 
: . ; 6 — 0.02264 : = 30 
0.1038 - - 373 cos 6 — 0.0226 ; — 378 cos 
57 ( : 30 0.03727 4 
0.06157 Br8 = 5rd cos + Udiad Sr8 —_ 5rs cos 
1 1 
36 + 0.003156 (: - +) cos 5@ + 0.03036 X 
r r 
( 1 1 ) ” 
—— = —- Joos. 
5r5 qr’ 
i iC “@ 10 d ( 1 
5. = — —- — M=ji1 - 3.8986 
its So: u r 00 (* j + r? + . 


“ 


1 1 1 1 ” q 
-_— — 4.8998 | — — =—] should read = 
3r4 5r6 5r6 7r8 u 


10/¢\, ‘) — ( 1 1 
- M= .4687 | —— 
r 2(*), od (1 + r? S r? = 3r4 - 


25 . : + 0.075 : z 4 
o.taee 3r4 5 r8 ee \ 58 778)" 
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Lines 31 and 32: “Ay,’ = 0.8307” should read “A;,’ = 


0.08307.” 


ooo, 


Induction Effects in Aeroelasticity—Errata 


John W. Miles 

Assistant Professor of Engineering, University of California at Los 
Angeles 

January, 1949 


(Readers’ Forum, 
1, pp. 
two figures were inadvertently omitted. 


In “Induction Effects in Aeroclasticity” 
JOURNAL OF THE AERONAUTICAL SCIENCES, Vol. 16, No 
63-64, January, 1949 
Figs. 1 and 2 are reproduced herewith, together with that part of 
the text which includes the reference to these figures as follows: 


In the case of control surface loading the reversal speed is more 
significant than the divergence speed, and it is therefore expedient 
to introduce the parameter @g as the ratio of the indicated air 
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speed to the reversal speed (@g = 1 at reversal). In this case it is 
found that the results are relatively insensitive to induction 
effects if the two-dimensional lift curve slope is used. The ratio 
(@p/@), as given by Eqs. (1) or (2) [with (dcz/da) evaluated as 
the two dimensional lift curve slope], is plotted against the 
parameter, I's in the accompanying Fig. 1 for a 30 per cent span 
aileron on a two to one tapered wing of aspect ratio 6, where 
TMs =1+ (cms /eé 15) (5) 
Cm, and C1, being the section characteristics of the control surface 
These results are not too sensitive to variations of parameters 
other than T's.- Using N.A.C.A. test data,’ Ig is plotted versus 
the control surface chord ratio for various products of the two 
dimensional lift curve slope (m) and the eccentricity (e) in Fig. 2 


* Ames, M. B., and Sears, R. I., Determination of Control Surface Char 


acteristics from N.A.C.A. Plain Flap and Tab Data, N.A.C.A. T.R. No. 721 
1941. 


On Chordwise Divergence* 


John W. Miles 

Assistant Professor of Engineering, University of California at Los 
Angeles 

November 22, 1948 


apd HAS RECENTLY CONSIDERED the problem of chordwise 
divergence of a two-dimensional airfoil at supersonic speeds 
on the assumption of a bending stiffness proportional to the third 
power of the thickness. If the airfoil is described in terms of the 
dimensionless coordinate £, running from @ at the leading edge 


to 1 at the section where it is assumed fixed (see Fig. 1), the bend 











| 
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Section of airfoil. 


Fic. 1. 


ing deflection due to the air loads satisfies the differential equa 


tion 
d? d*h(é) dh(£) 
T(é k— =0 I 
al (8) dz? | t dé 
k = 2M2(M? — 1)~"/* pc?/EI(1) 2 


where MM is the Mach Number, c is the velocity of sound, /(£) is 
the structurally effective moment of inertia relative to (1), and 
EI(£) is the bending stiffness per unit length. The problem is t 
find the eigenvalues of & such that Eq. (1) has a nonzero solution 
Biot considered the case where /(£) = £&, so that Eq. (1) is of the 
Euler type and yields solutions without resorting to infinite series 

Biot remarks that Eq. (1) is not self-adjoint, so that the familiar 
method of Rayleigh is not applicable (at least without further 
study and justification). A similar remark has been :nade by 
A. Flax? with respect to the spanwise bending of a swept wing 


under aerodynamic loading. 


* This work was carried out at North American Aviation, Inc. 

1 Biot, M. A., Aeroelastic Stability of Supersonic Wings—Chordwise Ds 
vergence, the Two-Dimensional Case, Cornell Aero. Lab. Report CM-427 
December 8, 1947. 

2 Sears, W. R., and Pai, S. I., Aeroelastic Effects for Swept Wings, 1. AS 
Preprint No 162, see footnote therein 
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rie. S k(1 —a)* vs.a for T(é) = £3, 


While it is true that the eigenvalues to the non-self-adjoint 
quation are no longer stationary, it does not follow that the 
Rayleigh method (and related methods) are no longer valid, and 
it is the writer’s opinion that the Galerkin form of Rayleigh’s 
method’ is quite suitable. To apply this method the bending 
jeflection is expanded in a set of functions y,,(£) satisfying canti- 


lever boundary conditions, viz., 
bni(1) = ¥,’(1) = T(a) van"(a) = [T(e)pn"(£)] ‘e ne (3) 


Substituting the expansion in Eq. (1), multiplying both sides of 
the equation by y,,(é)dé, integrating from — = a to — = 1, inte- 
grating the first term twice by parts, allowing m to assume all 
values assumed by , and requiring the determinant of the re- 
sulting set of homogeneous equations (in the expansion coeffi- 
cients) to vanish yields 
y ao | 

/ T(E) ym" (En (E)dzE + k Wm(£) Wn’ (t)dé = 0 (4) 
Ja 


a 


[fit is assumed that the set are the eigenfunctions, the eigenvalues 


“i 
2(a) / T(t) [Wn"(2)] 2d (5) 


In practice, a suitable, finite set of eigenfunctions may be esti- 


ire given by 


nated, and the divergence speed is specified by the smallest root 


*Duncan, W. J., The Principles of the Galerxin Method, R. & M. 1848, 


ARC, London, 1938. 
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k(1 — a)? vs. a for 7(£) = £2 
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In particular, if only one function is assumed, the 
It is physically evident 


to Eq. (4).‘ 
required eigenvalue is given by Eq. (5). 
that the true deflection function in the lowest mode lies some- 
where between the deflection functions for the tip loaded and 
uniformly loaded cantilever beam, and examination of Eq. (5) 
reveals that these functions will give lower and upper bounds, 
respectively, to the true eigenvalue. Evaluating these deflection 
functions for a cantilever beam of stiffness distribution 7(£) and 
substituting in Eq. (5) yields 


1 1 
k, =2 I-()( — a)*%dé (6) 
a 
"1 ~2 1 
ky = 2 I-'(&)(& — a)%dé I-(t)(£ — a)dé (7) 
a a 


for the lower and upper bounds. 

Comparing the results [Eqs. (6) and (7)] with Biot’s exact 
results! for the case /(£) = £ yields the curves shown in Fig. 2. 
In the case of a double wedge airfoil, a good approximation 
is T(t) = £2, and Eqs. (6) and (7) reduce to 


k, = 2[1 — a? + 2alna|— (8) 


l 
| (1 — a)(1 — 5a + 13a? + 3a) 4 tana | 
3 (9) 


ky = i : 
[; (1 — a) (1 — Sa — 2a’) - 2a%Ina 


These results are plotted in Fig. 3. 

Dr. C. C. Davenport (North American Aviation, Inc.) has ap- 
plied the Stodola (iteration) method to this problem. For T(t) = 
£3 his results compared with Biot’s after two iterations, and his 
results for 7(é) = £ are plotted in Fig. 3. Davenport has also 
considered the inhomogeneous problem, where the wing has an 
initial angle of attack, using the yield stress rather than diver- 
gence as a limiting criterion. 

The foregoing results are of interest both because of their direct 
ipplication and because of the indication of the accuracy to be 
expected in other non-self-adjoint problems, such as the bending 
of a swept wing, as treated by the writer. 

4 The Rayleigh principle implies that the convergence of the approximate 
eigenvalues on their true values is monotonic as the number of degrees of 
freedom is increased, by virtue of the associated variational principle. This 


result is not valid here. 
5 Miles, J. W., A Formulation of the Aeroelastic Problem f 
Northrop Aircraft Report GM-119 (A-96), May, 1948 


ra Swept Wing, 


Low-Speed Flutter 


J. H. Greidanus 


National Aeronautical Research Institute, Amsterdam, Holland 
November 29, 1948 


| gare ona TO THE INTERESTING PAPER ‘“‘Low-Speed Flutter 
and Its Physical Interpretation” (presented by M. A. Biot 
and Lee Arnold in the April, 1948, JouRNAL OF THE AERONAU- 
TICAL SCIENCES), it may be of interest to compare the conclu- 
sion, that the location of the nodal line is a dominating factor in 
the development of low-speed flutter, with the results of calcu- 
lations carried out a few years ago in the National Aeronautical 
Research Institute in Amsterdam, Holland, and published in a 
report! written in Dutch and probably not known outside my 
country. 

Considering the two-dimensional case of a rigid airfoil restrained 
by springs and capable of performing translational oscillations 
and rotational oscillations about the quarter-chord axis (not 
necessarily identical with the axis of application of the springs), 
we calculated the mean work £& per unit of time, done by the 

1! Greidanus, J. H., De arbeid, die de luchtkrachten per tijdseenheid op een 


trillende vleugel verrichten, ens., Report V. 1237 of the National Aeronautical 
Research Institute (NLL), 1940 
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aerodynamic forces under all possible kinematically different 


conditions of harmonic oscillation. We found that 


: l 
i = Tp 


1 


v3go*Cr 


1 
Ck Mya t)? + { fe(v) sin y + fs(v) cos WP} bo — 


where p = air density; ¢ = chord; v = frequency (rad. sec. ~'); 
v = reduced velocity = speed of flight divided by the product 
ve; & = dimensionless ratio y)/¢oc of the (real) amplitudes yo, 
¢o of translational and rotational kinematical components of the 
oscillation; y = phase lag of translational with respect to ro- 
tational component; fi, fe, f, = certain functions of v, based on 
the Kussner aerodynamic theory of the oscillating airfoil. 


All critical oscillations are, obviously, given by the equation 
Cr(v, fo, ¥) = O 
represented in Fig. 1. In this figure the area within the loops is 
unstable. Forv = ©, 1/2v = 0, the part0 < &,0< y¥ < x of 
the , ¥ plane—i.e., exactly one-half of the area 0 < h,0S y< 
2x covering all possible kinematically different combined oscilla- 
tions—becomes unstable. The importance of the amplitude 
ratio—i.e., of the location of the nodal line if the phase shift 
vanishes—is at once visible. The value & = 1/2, denoting ro- 
tation about the three-quarter-chord axis (y is taken positive 
upward and ¢ positive backward) if no phase shift occurs, is 
clearly the sole limiting case of actual flutter at zero speed of 
flight (and vanishing structural damping), and the same value of 
f) is, when combined with phase shift at any speed of flight (re- 
duced velocity v), unstable for a range of phase shifts of prac- 
tically maximum extent. There is, in the presence of phase 
shift, actually no nodal line but only an axis of minimum ampli- 
tude. Its distance xc aft of the quarter-chord line is determined 
by 
x = Scosy 

It coincides with the three-quarter-chord line if f) cos y = 1/2, a 
relation represented by a curve in Fig. 1, which also appears to 
pass right through the dangerous area. For a given value of », 
a certain ‘‘extremely unstable”’ oscillation can be defined by the 


Critical oscillations. 


maximum of the coefficient cg in the &, Y plane. (The case is, 
actually, an artificial one, since the calculation of EF does not 
apply to increasing oscillations but only to strictly harmoni¢ 
oscillations.) The pertaining point satisfies the equations 
Och /0E) = Oand Ocg/OY = 0. The locus of these points has also 
been introduced in the figure and is seen to pass somewhere be 
tween the & = 1/2 axis and the & cos y = 1/2 curve. The 
maximum value of cg, v given, increases linearly with v for sma 

values of v (up to about 0.6). It, however, soon increases mo 

quickly and shows a steep rise beyond v = 2 


Errata 


John W. Miles 

Assistant Professor of Engineering, University of California at Los 
Angeles 

November 19, 1948 


F E. Futter (Lockheed Aircraft Corporation) has brought 
e it to my attention that the formula for Cig due to fuselag 
in a recent paper! was incorrectly written because of a misinter 
pretation of Multhopp’s notation.? 

The corrected result should read: 


mo (2)(ie)(% + bp 
Cig = + ——— |“ ~—F)a 
1 + (2m cos o/7A.R.) 4 br] b 


(33) 


where the notation is that of reference 1, Eqs. (33), (34). Mr 
Fuller states that, in this form, the equation is in accord with 
known experimental data. 


1 Miles, J. W., The Rolling Moment Due to Sideslip for a Swept Wing 
Journal of the Aeronautical Sciences, Vol. 15, pp. 418-424, Eq. (33), July, 
1948. 

2? Multhopp, H., 
Eq. (7.10), December, 1942. 

3 Bamber, M. J., and House, R. O., Wind Tunnel Investigation of Effea 
of Yaw on Lateral Stability Characteristics, N.A.C.A. T.N. No. 730, 1939. 


Aerodynamics of the Fuselage, N.A.C.A. T.M. No. 1036, 
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